content: /*
   *    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., 675 Mass Ave, Cambridge, MA 02139, USA.
   */
  
  /*
   *    RegSMOImproved.java
   *    Copyright (C) 2006 University of Waikato, Hamilton, New Zealand
   *
   */
  
  package weka.classifiers.functions.supportVector;
  
  import weka.core.Instances;
  import weka.core.Option;
  import weka.core.RevisionUtils;
  import weka.core.TechnicalInformation;
  import weka.core.TechnicalInformationHandler;
  import weka.core.Utils;
  import weka.core.TechnicalInformation.Field;
  import weka.core.TechnicalInformation.Type;
  
  import java.util.Enumeration;
  import java.util.Vector;
  
  /**
   <!-- globalinfo-start -->
   * Learn SVM for regression using SMO with Shevade, Keerthi, et al. adaption of the stopping criterion.<br/>
   * <br/>
   * For more information see:<br/>
   * <br/>
   * S.K. Shevade, S.S. Keerthi, C. Bhattacharyya, K.R.K. Murthy: Improvements to the SMO Algorithm for SVM Regression. In: IEEE Transactions on Neural Networks, 1999.<br/>
   * <br/>
   * S.K. Shevade, S.S. Keerthi, C. Bhattacharyya, K.R.K. Murthy (1999). Improvements to the SMO Algorithm for SVM Regression. Control Division, Dept. of Mechanical Engineering.
   * <p/>
   <!-- globalinfo-end -->
   *
   <!-- technical-bibtex-start -->
   * BibTeX:
   * <pre>
   * &#64;inproceedings{Shevade1999,
   *    author = {S.K. Shevade and S.S. Keerthi and C. Bhattacharyya and K.R.K. Murthy},
   *    booktitle = {IEEE Transactions on Neural Networks},
   *    title = {Improvements to the SMO Algorithm for SVM Regression},
   *    year = {1999},
   *    PS = {http://guppy.mpe.nus.edu.sg/\~mpessk/svm/ieee_smo_reg.ps.gz}
   * }
   * 
   * &#64;techreport{Shevade1999,
   *    address = {Control Division, Dept. of Mechanical Engineering},
   *    author = {S.K. Shevade and S.S. Keerthi and C. Bhattacharyya and K.R.K. Murthy},
   *    institution = {National University of Singapore},
   *    number = {CD-99-16},
   *    title = {Improvements to the SMO Algorithm for SVM Regression},
   *    year = {1999},
   *    PS = {http://guppy.mpe.nus.edu.sg/\~mpessk/svm/smoreg_mod.ps.gz}
   * }
   * </pre>
   * <p/>
   <!-- technical-bibtex-end -->
   *
   <!-- options-start -->
   * Valid options are: <p/>
   * 
   * <pre> -T &lt;double&gt;
   *  The tolerance parameter for checking the stopping criterion.
   *  (default 0.001)</pre>
   * 
   * <pre> -V
   *  Use variant 1 of the algorithm when true, otherwise use variant 2.
   *  (default true)</pre>
   * 
   * <pre> -P &lt;double&gt;
   *  The epsilon for round-off error.
   *  (default 1.0e-12)</pre>
   * 
   * <pre> -L &lt;double&gt;
   *  The epsilon parameter in epsilon-insensitive loss function.
   *  (default 1.0e-3)</pre>
   * 
   * <pre> -W &lt;double&gt;
   *  The random number seed.
   *  (default 1)</pre>
   * 
   <!-- options-end -->
   *
   * @author  Remco Bouckaert (remco@cs.waikato.ac.nz,rrb@xm.co.nz)
   * @version $Revision: 1.4 $
   */
  public class RegSMOImproved
    extends RegSMO
    implements TechnicalInformationHandler {
    
    /** for serialization */
    private static final long serialVersionUID = 471692841446029784L;
    
    public final static int I0 = 3;
    public final static int I0a = 1;
    public final static int I0b = 2;
    public final static int I1 = 4;
    public final static int I2 = 8;
    public final static int I3 = 16;
    
    /** The different sets used by the algorithm. */
    protected SMOset m_I0;
    
    /** Index set {i: 0 < m_alpha[i] < C || 0 < m_alphaStar[i] < C}} */
    protected int [] m_iSet;
    
    /** b.up and b.low boundaries used to determine stopping criterion */
    protected double m_bUp, m_bLow;
    
    /** index of the instance that gave us b.up and b.low */
    protected int m_iUp, m_iLow;
    
    /** tolerance parameter used for checking stopping criterion b.up < b.low + 2 tol */
    double m_fTolerance = 0.001;
    
    /** set true to use variant 1 of the paper, otherwise use variant 2 */
    boolean m_bUseVariant1 = true;
    
    /**
     * Returns a string describing the object
     * 
     * @return 		a description suitable for
     * 			displaying in the explorer/experimenter gui
     */
    public String globalInfo() {
      return 
          "Learn SVM for regression using SMO with Shevade, Keerthi, et al. " 
        + "adaption of the stopping criterion.\n\n"
        + "For more information see:\n\n"
        + getTechnicalInformation().toString();
    }
  
    /**
     * Returns an instance of a TechnicalInformation object, containing 
     * detailed information about the technical background of this class,
     * e.g., paper reference or book this class is based on.
     * 
     * @return 		the technical information about this class
     */
    public TechnicalInformation getTechnicalInformation() {
      TechnicalInformation 	result;
      TechnicalInformation	additional;
      
      result = new TechnicalInformation(Type.INPROCEEDINGS);
      result.setValue(Field.AUTHOR, "S.K. Shevade and S.S. Keerthi and C. Bhattacharyya and K.R.K. Murthy");
      result.setValue(Field.TITLE, "Improvements to the SMO Algorithm for SVM Regression");
      result.setValue(Field.BOOKTITLE, "IEEE Transactions on Neural Networks");
      result.setValue(Field.YEAR, "1999");
      result.setValue(Field.PS, "http://guppy.mpe.nus.edu.sg/~mpessk/svm/ieee_smo_reg.ps.gz");
      
      additional = result.add(Type.TECHREPORT);
      additional.setValue(Field.AUTHOR, "S.K. Shevade and S.S. Keerthi and C. Bhattacharyya and K.R.K. Murthy");
      additional.setValue(Field.TITLE, "Improvements to the SMO Algorithm for SVM Regression");
      additional.setValue(Field.INSTITUTION, "National University of Singapore");
      additional.setValue(Field.ADDRESS, "Control Division, Dept. of Mechanical Engineering");
      additional.setValue(Field.NUMBER, "CD-99-16");
      additional.setValue(Field.YEAR, "1999");
      additional.setValue(Field.PS, "http://guppy.mpe.nus.edu.sg/~mpessk/svm/smoreg_mod.ps.gz");
      
      return result;
    }
    
    /**
     * Returns an enumeration describing the available options
     * 
     * @return 		an enumeration of all the available options
     */
    public Enumeration listOptions() {
      Vector result = new Vector();
      
      result.addElement(new Option(
  	"\tThe tolerance parameter for checking the stopping criterion.\n" 
  	+ "\t(default 0.001)", 
  	"T", 1, "-T <double>"));
      
      result.addElement(new Option(
  	"\tUse variant 1 of the algorithm when true, otherwise use variant 2.\n" 
  	+ "\t(default true)", 
  	"V", 0, "-V"));
      
      Enumeration enm = super.listOptions();
      while (enm.hasMoreElements()) {
        result.addElement(enm.nextElement());
      }
      
      return result.elements();
    }
    
    /**
     * Parses a given list of options. <p/>
     *
     <!-- options-start -->
     * Valid options are: <p/>
     * 
     * <pre> -T &lt;double&gt;
     *  The tolerance parameter for checking the stopping criterion.
     *  (default 0.001)</pre>
     * 
     * <pre> -V
     *  Use variant 1 of the algorithm when true, otherwise use variant 2.
     *  (default true)</pre>
     * 
     * <pre> -P &lt;double&gt;
     *  The epsilon for round-off error.
     *  (default 1.0e-12)</pre>
     * 
     * <pre> -L &lt;double&gt;
     *  The epsilon parameter in epsilon-insensitive loss function.
     *  (default 1.0e-3)</pre>
     * 
     * <pre> -W &lt;double&gt;
     *  The random number seed.
     *  (default 1)</pre>
     * 
     <!-- options-end -->
     *
     * @param options 	the list of options as an array of strings
     * @throws Exception 	if an option is not supported 
     */
    public void setOptions(String[] options) throws Exception {
      String	tmpStr;
      
      tmpStr = Utils.getOption('T', options);
      if (tmpStr.length() != 0) {
        setTolerance(Double.parseDouble(tmpStr));
      } else {
        setTolerance(0.001);
      }
      
      setUseVariant1(Utils.getFlag('V', options));
      
      super.setOptions(options);
    }
    
    /**
     * Gets the current settings of the object.
     *
     * @return 		an array of strings suitable for passing to setOptions
     */
    public String[] getOptions() {
      int       	i;
      Vector    	result;
      String[]  	options;
  
      result = new Vector();
  
      options = super.getOptions();
      for (i = 0; i < options.length; i++)
        result.add(options[i]);
      
      result.add("-T");
      result.add("" + getTolerance());
      
      if (m_bUseVariant1)
        result.add("-V");
  
      return (String[]) result.toArray(new String[result.size()]);	  
    }
    
    /**
     * Returns the tip text for this property
     * 
     * @return 		a description suitable for
     * 			displaying in the explorer/experimenter gui
     */
    public String toleranceTipText() {
      return "tolerance parameter used for checking stopping criterion b.up < b.low + 2 tol";
    }
    
    /**
     * returns the current tolerance
     * 
     * @return	the tolerance
     */
    public double getTolerance() {
      return m_fTolerance;
    }
    
    /**
     * sets the tolerance
     * 
     * @param d	the new tolerance
     */
    public void setTolerance(double d) {
      m_fTolerance = d;
    }
    
    /**
     * Returns the tip text for this property
     * 
     * @return 		a description suitable for
     * 			displaying in the explorer/experimenter gui
     */
    public String useVariant1TipText() {
      return "set true to use variant 1 of the paper, otherwise use variant 2.";
    }
    
    /**
     * Whether variant 1 is used
     * 
     * @return		true if variant 1 is used
     */
    public boolean isUseVariant1() {
      return m_bUseVariant1;
    }
    
    /**
     * Sets whether to use variant 1
     * 
     * @param b		if true then variant 1 is used
     */
    public void setUseVariant1(boolean b) {
      m_bUseVariant1 = b;
    }
    
    /** 
     * takeStep method from Shevade et al.s paper.
     * parameters correspond to pseudocode from paper.
     * 
     * @param i1
     * @param i2
     * @param alpha2
     * @param alpha2Star
     * @param phi2
     * @return
     * @throws Exception
     */
    protected int takeStep(int i1, int i2, double alpha2, double alpha2Star, double phi2) throws Exception {
      //procedure takeStep(i1, i2)
      //
      //  if (i1 == i2) 
      //    return 0 
      if (i1 == i2) {
        return 0;
      }
      double C1 = m_C * m_data.instance(i1).weight();
      double C2 = m_C * m_data.instance(i2).weight();
      //  alpha1, alpha1' = Lagrange multipliers for i1 
      double alpha1 = m_alpha[i1];
      double alpha1Star = m_alphaStar[i1];
  //  double y1 = m_target[i1];
      // TODO: verify we do not need to recompute m_error[i1] here
      // TODO: since m_error is only updated for indices in m_I0
      double phi1 = m_error[i1];
  //  if ((m_iSet[i1] & I0)==0) {
  //  phi1 = -SVMOutput(i1) - m_b + m_target[i1];
  //  m_error[i1] = phi1;
  //  }
      //  k11 = kernel(point[i1], point[i1]) 
      //  k12 = kernel(point[i1], point[i2]) 
      //  k22 = kernel(point[i2], point[i2]) 
      //  eta = -2*k12+k11+k22 
      //  gamma = alpha1-alpha1'+alpha2-alpha2'
      //
      double k11 = m_kernel.eval(i1, i1, m_data.instance(i1));
      double k12 = m_kernel.eval(i1, i2, m_data.instance(i1));
      double k22 = m_kernel.eval(i2, i2, m_data.instance(i2));
      double eta = -2 * k12 + k11 + k22;
      double gamma = alpha1 - alpha1Star + alpha2 - alpha2Star;
  //  if (eta < 0) {
      // this may happen due to numeric instability
      // due to Mercer's condition, this should not happen, hence we give up
  //  return 0;
  //  }
      //  % We assume that eta > 0. Otherwise one has to repeat the complete 
      //  % reasoning similarly (i.e. compute objective functions at L and H 
      //  % and decide which one is largest
      //
      //  case1 = case2 = case3 = case4 = finished = 0 
      //  alpha1old = alpha1, 
      //  alpha1old' = alpha1' 
      //  alpha2old = alpha2, 
      //  alpha2old' = alpha2' 
      //  deltaphi = F1 - F2 
      //
      
      //  while !finished
      //    % This loop is passed at most three times 
      //    % Case variables needed to avoid attempting small changes twice 
      //    if (case1 == 0) &&
      //       (alpha1 > 0 || (alpha1' == 0 && deltaphi > 0)) && 
      //       (alpha2 > 0 || (alpha2' == 0 && deltaphi < 0))
      //        compute L, H (w.r.t. alpha1, alpha2) 
      //        if (L < H)
      //          a2 = alpha2 - (deltaphi / eta ) a2 = min(a2, H) a2 = max(L, a2) a1 = alpha1 - (a2 - alpha2) 
      //          update alpha1, alpha2 if change is larger than some eps 
      //        else
      //          finished = 1 
      //        endif 
      //      case1 = 1 
      //    elseif (case2 == 0) &&
      //           (alpha1 > 0 || (alpha1' == 0 && deltaphi > 2*epsilon)) && 
      //           (alpha2' > 0 || (alpha2 == 0 && deltaphi > 2*epsilon))
      //
      //        compute L, H (w.r.t. alpha1, alpha2') 
      //        if (L < H)
      //          a2 = alpha2' + ((deltaphi - 2*epsilon)/eta)) a2 = min(a2, H) a2 = max(L, a2) a1 = alpha1 + (a2-alpha2') 
      //          update alpha1, alpha2' if change is larger than some eps 
      //        else
      //          finished = 1 
      //        endif 
      //        case2 = 1 
      //    elseif (case3 == 0) &&
      //           (alpha1' > 0 || (alpha1 == 0 && deltaphi < -2*epsilon)) && 
      //           (alpha2 > 0 || (alpha2' == 0 && deltaphi < -2*epsilon))
      //         compute L, H (w.r.t. alpha1', alpha2) 
      //         if (L < H)
      //           a2 = alpha2 - ((deltaphi + 2*epsilon)/eta) a2 = min(a2, H) a2 = max(L, a2) a1 = alpha1' + (a2 - alpha2) 
      //           update alpha1', alpha2 if change is larger than some eps 
      //         else
      //           finished = 1 
      //         endif 
      //         case3 = 1 
      //    elseif (case4 == 0) &&
      //           (alpha1' > 0) || (alpha1 == 0 && deltaphi < 0)) && 
      //           (alpha2' > 0) || (alpha2 == 0 && deltaphi > 0))
      //         compute L, H (w.r.t. alpha1', alpha2') 
      //         if (L < H) 
      //           a2 = alpha2' + deltaphi/eta a2 = min(a2, H) a2 = max(L, a2) a1 = alpha1' - (a2 - alpha2') 
      //           update alpha1, alpha2' if change is larger than some eps 
      //         else
      //           finished = 1 
      //         endif 
      //         case4 = 1 
      //    else
      //      finished = 1 
      //    endif 
      //    update deltaphi 
      //  endwhile 
      
      double alpha1old = alpha1;
      double alpha1Starold = alpha1Star;
      double alpha2old = alpha2;
      double alpha2Starold = alpha2Star;
      double deltaPhi = phi1 - phi2;
      
      if (findOptimalPointOnLine(i1, alpha1, alpha1Star, C1, i2, alpha2, alpha2Star, C2, gamma, eta, deltaPhi)) {
        
        alpha1 = m_alpha[i1];
        alpha1Star = m_alphaStar[i1];
        alpha2 = m_alpha[i2];
        alpha2Star = m_alphaStar[i2];
        
        //  if changes in alpha('), alpha2(') are larger than some eps
        //    Update f-cache[i] for i in I.0 using new Lagrange multipliers 
        //    Store the changes in alpha, alpha' array 
        //    Update I.0, I.1, I.2, I.3 
        //    Compute (i.low, b.low) and (i.up, b.up) by applying the conditions mentioned above, using only i1, i2 and indices in I.0 
        //    return 1 
        //  else
        //    return 0
        //endif endprocedure
        
        //		Update error cache using new Lagrange multipliers 
        double dAlpha1 = alpha1 - alpha1old - (alpha1Star - alpha1Starold);
        double dAlpha2 = alpha2 - alpha2old - (alpha2Star - alpha2Starold);
        for (int j = m_I0.getNext(-1); j != -1; j = m_I0.getNext(j)) {
  	if ((j != i1) && (j != i2)) {
  	  m_error[j] -= dAlpha1 * m_kernel.eval(i1, j, m_data.instance(i1)) 
  	  + dAlpha2 * m_kernel.eval(i2, j, m_data.instance(i2));
  	}
        }
        m_error[i1] -= dAlpha1 * k11 + dAlpha2 * k12;
        m_error[i2] -= dAlpha1 * k12 + dAlpha2 * k22;
        
        updateIndexSetFor(i1, C1);
        updateIndexSetFor(i2, C2);
        
        //    Compute (i.low, b.low) and (i.up, b.up) by applying the conditions mentioned above, using only i1, i2 and indices in I.0 
        m_bUp = Double.MAX_VALUE; 
        m_bLow = -Double.MAX_VALUE; 
        for (int j = m_I0.getNext(-1); j != -1; j = m_I0.getNext(j)) {
  	updateBoundaries(j, m_error[j]);
        }
        if (!m_I0.contains(i1)) {
  	updateBoundaries(i1, m_error[i1]);
        }
        if (!m_I0.contains(i2)) {
  	updateBoundaries(i2, m_error[i2]);
        }
        
        return 1;
      } 
      else {
        return 0;
      }
    }
    
    /**
     * updates the index sets I0a, IOb, I1, I2 and I3 for vector i
     * 
     * @param i index of vector
     * @param C capacity for vector i
     * @throws Exception
     */
    protected void updateIndexSetFor(int i, double C) throws Exception {
      /*
       m_I0a.delete(i);
       m_I0b.delete(i);
       m_I1.delete(i);
       m_I2.delete(i);
       m_I3.delete(i);
       */
      if (m_alpha[i] == 0 && m_alphaStar[i] == 0) {
        //m_I1.insert(i);
        m_iSet[i] = I1;
        m_I0.delete(i);
      } else if (m_alpha[i] > 0) {
        if (m_alpha[i] < C) {
  	if ((m_iSet[i] & I0) == 0) {
  	  //m_error[i] = -SVMOutput(i) - m_b + m_target[i];
  	  m_I0.insert(i);
  	}
  	//m_I0a.insert(i);
  	m_iSet[i] = I0a;
        } else { // m_alpha[i] == C
  	//m_I3.insert(i);
  	m_iSet[i] = I3;
  	m_I0.delete(i);
        }
      } else {// m_alphaStar[i] > 0 
        if (m_alphaStar[i] < C) {
  	if ((m_iSet[i] & I0) == 0) {
  	  //m_error[i] = -SVMOutput(i) - m_b + m_target[i];
  	  m_I0.insert(i);
  	}
  	//m_I0b.insert(i);
  	m_iSet[i] = I0b;
        } else { // m_alpha[i] == C
  	//m_I2.insert(i);
  	m_iSet[i] = I2;
  	m_I0.delete(i);
        }
      }
    }
    
    /**
     * updates boundaries bLow and bHi and corresponding indexes
     * 
     * @param i2 index of vector
     * @param F2 error of vector i2
     */
    protected void updateBoundaries(int i2, double F2) {		
      int iSet = m_iSet[i2];
      
      double FLow = m_bLow;
      if ((iSet & (I2 | I0b)) > 0) {
        FLow = F2 + m_epsilon;
      } else if ((iSet & (I1 | I0a)) > 0) {
        FLow = F2 - m_epsilon;
      }
      if (m_bLow < FLow) {
        m_bLow = FLow;
        m_iLow = i2;
      }
      double FUp = m_bUp;
      if ((iSet & (I3 | I0a)) > 0) {
        FUp = F2 - m_epsilon;
      } else if ((iSet & (I1 | I0b)) > 0) {
        FUp = F2 + m_epsilon;
      }
      if (m_bUp > FUp) {
        m_bUp = FUp;
        m_iUp = i2;
      }
    }
    
    /** 
     * parameters correspond to pseudocode from paper.
     * 
     * @param i2 index of  candidate
     * @return
     * @throws Exception
     */
    protected int examineExample(int i2) throws Exception {
      //procedure examineExample(i2)
      //
      //  alpha2, alpha2' = Lagrange multipliers for i2 
      double alpha2 = m_alpha[i2];
      double alpha2Star = m_alphaStar[i2];
      
      //  if (i2 is in I.0)
      //    F2 = f-cache[i2] 
      //  else
      //    compute F2 = F.i2 and set f-cache[i2] = F2 
      //    % Update (b.low, i.low) or (b.up, i.up) using (F2, i2)... 
      //    if (i2 is in I.1)
      //      if (F2+epsilon < b.up)
      //        b.up = F2+epsilon, 
      //        i.up = i2 
      //      elseif (F2-epsilon > b.low)
      //        b.low = F2-epsilon, 
      //        i.low = i2 
      //      end if 
      //    elseif ( (i2 is in I.2) && (F2+epsilon > b.low) )
      //      b.low = F2+epsilon, 
      //      i.low = i2 
      //    elseif ( (i2 is in I.3) && (F2-epsilon < b.up) )
      //      b.up = F2-epsilon, 
      //      i.up = i2 
      //    endif 
      //  endif 
      
      int iSet = m_iSet[i2];
      double F2 = m_error[i2];
      if (!m_I0.contains(i2)) {
        F2 = -SVMOutput(i2) - m_b + m_target[i2];
        m_error[i2] = F2;
        if (iSet == I1) {
  	if (F2 + m_epsilon < m_bUp) {
  	  m_bUp = F2 + m_epsilon;
  	  m_iUp = i2;
  	} else if (F2 - m_epsilon > m_bLow) {
  	  m_bLow = F2 - m_epsilon;
  	  m_iLow = i2;
  	} 
        } else if ((iSet == I2) && (F2 + m_epsilon > m_bLow)) {
  	m_bLow = F2 + m_epsilon;
  	m_iLow = i2;
        } else if ((iSet == I3) && (F2 - m_epsilon < m_bUp)) {
  	m_bUp = F2 - m_epsilon;
  	m_iUp = i2;
        }
      }
      
      //  % Check optimality using current b.low and b.up and, if 
      //  % violated, find an index i1 to do joint optimization with i2... 
      //  optimality = 1;
      //  case 1: i2 is in I.0a
      //    if (b.low-(F2-epsilon) > 2 * tol)
      //      optimality = 0;  
      //      i1 = i.low; 
      //      % For i2 in I.0a choose the better i1... 
      //      if ((F2-epsilon)-b.up > b.low-(F2-epsilon))
      //        i1 = i.up; 
      //      endif 
      //    elseif ((F2-epsilon)-b.up > 2 * tol)
      //      optimality = 0; 
      //      i1 = i.up; 
      //      % For i2 in I.0a choose the better i1... 
      //      if ((b.low-(F2-epsilon) > (F2-epsilon)-b.up)
      //        i1 = i.low; 
      //      endif 
      //    endif 
      //  case 2: i2 is in I.0b
      //    if (b.low-(F2+epsilon) > 2 * tol)
      //      optimality = 0; 
      //      i1 = i.low; 
      //      % For i2 in I.0b choose the better i1... 
      //      if ((F2+epsilon)-b.up > b.low-(F2+epsilon))
      //        i1 = i.up; 
      //      endif 
      //    elseif ((F2+epsilon)-b.up > 2 * tol)
      //      optimality = 0; 
      //      i1 = i.up; 
      //      % For i2 in I.0b choose the better i1... 
      //      if ((b.low-(F2+epsilon) > (F2+epsilon)-b.up)
      //        i1 = i.low; 
      //      endif 
      //    endif 
      //  case 3: i2 is in I.1
      //    if (b.low-(F2+epsilon) > 2 * tol)
      //      optimality = 0; 
      //      i1 = i.low; 
      //      % For i2 in I1 choose the better i1... 
      //      if ((F2+epsilon)-b.up > b.low-(F2+epsilon)
      //        i1 = i.up; 
      //      endif 
      //    elseif ((F2-epsilon)-b.up > 2 * tol)
      //      optimality = 0; 
      //      i1 = i.up; 
      //      % For i2 in I1 choose the better i1... 
      //      if (b.low-(F2-epsilon) > (F2-epsilon)-b.up)
      //        i1 = i.low; 
      //      endif 
      //    endif 
      //  case 4: i2 is in I.2
      //    if ((F2+epsilon)-b.up > 2*tol)
      //      optimality = 0, 
      //      i1 = i.up 
      //     endif 
      //  case 5: i2 is in I.3
      //    if ((b.low-(F2-epsilon) > 2*tol)
      //      optimality = 0, i1 = i.low 
      //    endif
      
      int i1 = i2;
      boolean bOptimality = true;
      //case 1: i2 is in I.0a
      if (iSet == I0a) {
        if (m_bLow - (F2 - m_epsilon) > 2 * m_fTolerance) {
  	bOptimality = false;
  	i1 = m_iLow;
  	//% For i2 in I .0 a choose the better i1...
  	if ((F2 - m_epsilon) - m_bUp > m_bLow - (F2 - m_epsilon)) {
  	  i1 = m_iUp;
  	}
        } else if ((F2 - m_epsilon) - m_bUp > 2 * m_fTolerance) {
  	bOptimality = false;
  	i1 = m_iUp;
  	//% For i2 in I.0a choose the better i1... 
  	if (m_bLow - (F2 - m_epsilon) > (F2 - m_epsilon) - m_bUp) {
  	  i1 = m_iLow;
  	}
        }
      } // case 2: i2 is in I.0b
      else if (iSet == I0b) {
        if (m_bLow - (F2 + m_epsilon) > 2 * m_fTolerance) {
  	bOptimality = false;
  	i1 = m_iLow; // % For i2 in I.0b choose the better i1... 
  	if ((F2 + m_epsilon) - m_bUp > m_bLow - (F2 + m_epsilon)) {
  	  i1 = m_iUp;
  	}
        } else if ((F2 + m_epsilon) - m_bUp > 2 * m_fTolerance) {
  	bOptimality = false;
  	i1 = m_iUp; // % For i2 in I.0b choose the better i1... 
  	if (m_bLow - (F2 + m_epsilon) > (F2 + m_epsilon) - m_bUp) {
  	  i1 = m_iLow;
  	}
        }
      } // case 3: i2 is in I.1
      else if (iSet == I1) {
        if (m_bLow - (F2 + m_epsilon) > 2 * m_fTolerance) {
  	bOptimality = false;
  	i1 = m_iLow;
  	//% For i2 in I1 choose the better i1...
  	if ((F2 + m_epsilon) - m_bUp > m_bLow - (F2 + m_epsilon)) {
  	  i1 = m_iUp;
  	}
        } else if ((F2 - m_epsilon) - m_bUp > 2 * m_fTolerance) {
  	bOptimality = false;
  	i1 = m_iUp; // % For i2 in I1 choose the better i1... 
  	if (m_bLow - (F2 - m_epsilon) > (F2 - m_epsilon) - m_bUp) {
  	  i1 = m_iLow;
  	}
        }
      } //case 4: i2 is in I.2
      else if (iSet == I2) {
        if ((F2 + m_epsilon) - m_bUp > 2 * m_fTolerance) {
  	bOptimality = false;
  	i1 = m_iUp;
        }
      } //case 5: i2 is in I.3
      else if (iSet == I3) {
        if (m_bLow - (F2 - m_epsilon) > 2 * m_fTolerance) {
  	bOptimality = false;
  	i1 = m_iLow;
        }
      }
      // if (optimality == 1) 
      //    return 0
      //  if (takeStep(i1, i2))
      //    return 1 
      //   else
      //    return 0 
      //  endif 
      //endprocedure
      if (bOptimality) {
        return 0;
      }
      return takeStep(i1, i2, m_alpha[i2], m_alphaStar[i2], F2);
    }
    
    /** 
     * initialize various variables before starting the actual optimizer 
     * 
     * @param data 	data set used for learning
     * @throws Exception	if something goes wrong
     */
    protected void init(Instances data) throws Exception {
      super.init(data);
      // from Keerthi's pseudo code:
      //  set alpha and alpha' to zero for every example set I.1 to contain all the examples 
      //  Choose any example i from the training set. 
      //  set b.up = target[i]+epsilon 
      //  set b.low = target[i]-espilon 
      //  i.up = i.low = i; 
      // Initialize sets
      m_I0 = new SMOset(m_data.numInstances());
      m_iSet = new int [m_data.numInstances()];
      for (int i = 0; i < m_nInstances; i++) {
        m_iSet[i] = I1;
      }
      // m_iUp = m_random.nextInt(m_nInstances);
      m_iUp = 0;
      m_bUp = m_target[m_iUp] + m_epsilon;
      m_iLow = m_iUp;
      m_bLow = m_target[m_iLow] - m_epsilon;
      //init error cache
      m_error = new double[m_nInstances];
      for (int i = 0; i < m_nInstances; i++) {
        m_error[i] = m_target[i];
      }
    }
    
    /** 
     * use variant 1 of Shevade's et al.s paper
     * 
     * @throws Exception	if something goes wrong
     */
    protected void optimize1() throws Exception {
      //% main routine for modification 1 procedure main
      //  while (numChanged > 0 || examineAll)
      //    numChanged = 0; 
      int nNumChanged = 0;
      boolean bExamineAll = true;
      //  while (numChanged > 0 || examineAll)
      //    numChanged = 0; 
      while (nNumChanged > 0 || bExamineAll) {
        nNumChanged = 0;
        //    if (examineAll)
        //      loop I over all the training examples
        //        numChanged += examineExample(I) 
        //    else
        //      loop I over I.0
        //        numChanged += examineExample(I) 
        //        % It is easy to check if optimality on I.0 is attained... 
        //        if (b.up > b.low - 2*tol) at any I
        //          exit the loop after setting numChanged = 0 
        //        endif 
        if (bExamineAll) {
  	for (int i = 0; i < m_nInstances; i++) {
  	  nNumChanged += examineExample(i);
  	}
        } else {
  	for (int i = m_I0.getNext(-1); i != -1; i = m_I0.getNext(i)) {
  	  
  	  nNumChanged += examineExample(i);
  	  if (m_bLow - m_bUp < 2 * m_fTolerance) {
  	    nNumChanged = 0;
  	    break;
  	  }
  	}
        } //    if (examineAll == 1)
        //      examineAll = 0; 
        //    elseif (numChanged == 0)
        //      examineAll = 1;
        //    endif 
        //  endwhile 
        //endprocedure
        if (bExamineAll) {
  	bExamineAll = false;
        } else if (nNumChanged == 0) {
  	bExamineAll = true;
        }
      }
    }
    
    /** 
     * use variant 2 of Shevade's et al.s paper 
     * 
     * @throws Exception	if something goes wrong
     */
    protected void optimize2() throws Exception {
      //% main routine for modification 2 procedure main
      int nNumChanged = 0;
      boolean bExamineAll = true;
      //  while (numChanged > 0 || examineAll)
      //    numChanged = 0; 
      while (nNumChanged > 0 || bExamineAll) {
        nNumChanged = 0;
        //    if (examineAll)
        //      loop I over all the training examples
        //        numChanged += examineExample(I) 
        //    else
        //      % The following loop is the only difference between the two 
        //      % SMO modifications. Whereas, modification 1, the type II 
        //      % loop selects i2 fro I.0 sequentially, here i2 is always 
        //      % set to the current i.low and i1 is set to the current i.up; 
        //      % clearly, this corresponds to choosing the worst violating 
        //      % pair using members of I.0 and some other indices
        //      inner.loop.success = 1; 
        //      do
        //        i2 = i.low 
        //	      alpha2, alpha2' = Lagrange multipliers for i2 
        //	      F2 = f-cache[i2] 
        //	      i1 = i.up 
        //	      inner.loop.success = takeStep(i.up, i.low) 
        //	      numChanged += inner.loop.success 
        //      until ( (b.up > b.low - 2*tol) || inner.loop.success == 0) 
        //      numChanged = 0; 
        //    endif 
        if (bExamineAll) {
  	for (int i = 0; i < m_nInstances; i++) {
  	  nNumChanged += examineExample(i);
  	}
        } else {
  	boolean bInnerLoopSuccess = true;
  	do {
  	  if (takeStep(m_iUp, m_iLow, m_alpha[m_iLow], m_alphaStar[m_iLow], m_error[m_iLow]) > 0) {
  	    bInnerLoopSuccess = true;
  	    nNumChanged += 1;
  	  } else {
  	    bInnerLoopSuccess = false;
  	  }
  	} while ((m_bUp <= m_bLow - 2 * m_fTolerance) && bInnerLoopSuccess);
  	nNumChanged = 0;
        } //
        //    if (examineAll == 1)
        //      examineAll = 0 
        //    elseif (numChanged == 0)
        //      examineAll = 1 
        //    endif 
        //  endwhile 
        //endprocedure
        //
        if (bExamineAll) {
  	bExamineAll = false;
        } else if (nNumChanged == 0) {
  	bExamineAll = true;
        }
      }
    }
    
    /** 
     * wrap up various variables to save memeory and do some housekeeping after optimization
     * has finished.
     *
     * @throws Exception 	if something goes wrong
     */
    protected void wrapUp() throws Exception {
      m_b = -(m_bLow + m_bUp) / 2.0;
      m_target = null;
      m_error = null;
      super.wrapUp();
    }
    
    /** 
     * learn SVM parameters from data using Keerthi's SMO algorithm.
     * Subclasses should implement something more interesting.
     * 
     * @param instances	the data to work with
     * @throws Exception	if something goes wrong
     */
    public void buildClassifier(Instances instances) throws Exception {
      // initialize variables		
      init(instances); 
  
      // solve optimization problem
      if (m_bUseVariant1) {
        optimize1();
      } else {
        optimize2();
      } 
      
      // clean up
      wrapUp();
    }
    
    /**
     * Returns the revision string.
     * 
     * @return		the revision
     */
    public String getRevision() {
      return RevisionUtils.extract("$Revision: 1.4 $");
    }
  }
  