# Copyright 2000 by Jeffrey Chang.  All rights reserved.
# This code is part of the Biopython distribution and governed by its
# license.  Please see the LICENSE file that should have been included
# as part of this package.

"""This provides code for working with support vector machines.

Support vector machine (SVM) is a supervised machine learning
algorithm that has been shown to work well on a variety of
classification problems.

XXX describe kernel_fn, C, epsilon

For more information, see:
http://svm.first.gmd.de/
http://svm.research.bell-labs.com/

Classes:
SVM                  Holds data for a Support Vector Machine.
SMOTrainer           Trains an SVM using Sequential Minimal Optimization (SMO).
KeerthiTrainer       Trains an SVM using Keerthi's extensions to SMO.
TransductiveTrainer  Trains a transductive SVM.

LinearKernel
PolynomialKernel
RadialBasisFunctionKernel
HyperbolicTangentKernel


Functions:
train                Train a Support Vector Machine on some training data.
trans_transductive
classify             Use a Support Vector Machine to classify some data.

Usage:
The 'train' function is provided as a user-friendly interface module.
Use this function if all you want is a plain-vanilla SVM classifier.
However, if you're concerned about details such as the training
algorithm, you will need to use some of the classes provided.

"""
import warnings
warnings.warn("Bio.SVM is deprecated. We recommend using LIBSVM: "
              "http://www.csie.ntu.edu.tw/~cjlin/libsvm/",
              DeprecationWarning)

# To do:
# make update_fn better
#    take an event string, followed by data
#    or, make it an object that responds to various things

import math
import random

class LinearKernel:
    """k(x, y) = x*y"""
    def __init__(self):
        pass
    def __call__(self, x1, x2):
        return _dot(x1, x2)

# XXX document sparse kernels
# XXX try to generalize the kernels, so they work with either type of input
class SparseLinearKernel:
    """k(x, y) = x*y"""
    def __init__(self):
        pass
    def __call__(self, x1, x2):
        return _sparse_dot(x1, x2)

class PolynomialKernel:
    """k(x, y) = (x*y+1)**p"""
    _max_p = 20
    def __init__(self, p=3):
        if p < 1 or p > self._max_p:
            raise ValueError, \
                  "polynomial degree must be between 1 and %d" % self._max_p
        self.p = p
    def __call__(self, x1, x2):
        s = _dot(x1, x2) + 1.0
        return s ** self.p

class SparsePolynomialKernel:
    """k(x, y) = (x*y+1)**p"""
    _max_p = 20
    def __init__(self, p=3):
        if p < 1 or p > self._max_p:
            raise ValueError, \
                  "polynomial degree must be between 1 and %d" % self._max_p
        self.p = p
    def __call__(self, x1, x2):
        s = _sparse_dot(x1, x2) + 1.0
        return s ** self.p

class RadialBasisFunctionKernel:
    """      -(x-y)**2
         e** ----------
             2*sigma**2
    """
    def __init__(self, sigma=1.0):
        self.sigma = sigma
        self.sigma2 = sigma**2
    def __call__(self, x1, x2):
        x1_x2 = _subtract(x1, x2)
        num = -_dot(x1_x2, x1_x2)
        den = 2.0 * self.sigma2
        try:
            return math.exp(num/den)
        except OverflowError:
            return 0.0

class SparseRadialBasisFunctionKernel:
    """      -(x-y)**2
         e** ----------
             2*sigma**2
    """
    def __init__(self, sigma=1.0):
        self.sigma = sigma
        self.sigma2 = sigma**2
    def __call__(self, x1, x2):
        x1_x2 = _sparse_subtract(x1, x2)
        num = -_sparse_dot(x1_x2, x1_x2)
        den = 2.0 * self.sigma2
        try:
            return math.exp(num/den)
        except OverflowError:
            return 0.0

class HyperbolicTangentKernel:
    """tanh(x*y*kappa - delta)"""
    # not all kappa and delta obey Mercer's condition
    def __init__(self, kappa, delta):
        self.kappa = kappa
        self.delta = delta
    def __call__(self, x1, x2):
        x1x2 = _dot(x1, x2)
        return math.tanh(x1x2*self.kappa - self.delta)

class SVM:
    """Holds information for a non-linear Support Vector Machine.

    Members:
    xs          A list of the input vectors.
    ys          A list of the output values.  Must be either 1 or -1.
    alphas      A list of the Lagrange multipliers for each point.
    b           The threshold value for the machine.
    kernel_fn   Should take 2 vectors and return a distance.

    * xs, ys, and alphas should be a parallel list of vectors.
    ** There's a special read-only member variable 'w'.  For linear
    SVM's, this is the 'w' vector so that the decision hyperplane is at
    wx-b=0

    """
    def __init__(self, xs, ys, alphas=None, b=None, kernel_fn=LinearKernel()):
        if alphas is None:
            alphas = [0.0]*len(self.xs)
        if b is None:
            b = 0.0

        self.xs, self.ys = xs, ys
        self.alphas = alphas
        self.b = b
        self.kernel_fn = kernel_fn

        if len(xs) != len(ys):
            raise ValueError, "xs and ys must be parallel arrays"
        if len(xs) != len(alphas):
            raise ValueError, "xs and alphas must be parallel arrays"

        for i in range(len(alphas)):
            if self.ys[i] not in [1, -1]:
                raise ValueError, "Output must be either 1 or -1"

    def __getattr__(self, x):
        if x == 'w':
            w = [0.0] * len(self.xs[0])
            for i in range(len(w)):
                for j in range(len(self.xs)):
                    w[i] = w[i] + self.ys[j]*self.alphas[j]*self.xs[j][i]
            return w
        raise AttributeError, x
        
def classify(svm, x):
    """classify(svm, x) -> score

    Classify x based on an SVM object.

    """
    sum = 0.0
    for i in range(len(svm.xs)):
        if svm.alphas[i] == 0.0:
            continue
        sum = sum + svm.ys[i]*svm.alphas[i]*svm.kernel_fn(svm.xs[i], x)
    return sum - svm.b

class SMOTrainer:
    """Sequential Minimal Optimization Trainer

    This is an implementation of the SVM trainer called Sequential
    Minimal Optimization, described by John C. Platt:
    http://www.research.microsoft.com/~jplatt

    Methods:
    train     Train a new SVM.

    """
    def train(self, training_set, results, kernel_fn, C, epsilon,
              update_fn=None):
        """train(self, training_set, results, kernel_fn, C, epsilon,
        update_fn=None) -> SVM

        """
        # Initializing variables
        self.svm = SVM(training_set, results,
                       [0.0]*len(training_set), 0.0, kernel_fn)
        self.eps = epsilon
        self.C = C
        self._ecache = {}

        # Training algorithm starts here.
        num_changed = 0
        examine_all = 1

        while num_changed > 0 or examine_all:
            num_changed = 0
            if examine_all:
                for i2 in range(len(self.svm.alphas)):
                    num_changed = num_changed + self._examine_example(i2)
            else:
                for i2 in range(len(self.svm.alphas)):
                    # skip the bound alphas
                    if self._is_bound(self.svm.alphas[i2], self.C):
                        continue
                    num_changed = num_changed + self._examine_example(i2)
            if update_fn is not None:
                update_fn(num_changed, self.svm)
            if examine_all:
                examine_all = 0
            elif num_changed == 0:
                examine_all = 1

        return self.svm

    def _examine_example(self, i2):
        y2 = self.svm.ys[i2]
        alph2 = self.svm.alphas[i2]
        E2 = self._E(i2)
        r2 = E2*y2

        if (r2 < -self.eps and alph2 < self.C) or \
           (r2 > self.eps and alph2 > 0):
            # sort it by E1 - E2
            max_i1 = None
            max_error = None
            for i1 in range(len(self.svm.alphas)):
                if i1 == i2 or self._is_bound(self.svm.alphas[i1], self.C):
                    continue
                E1 = self._E(i1)
                error = math.fabs(E1 - E2)
                if max_i1 is None or error > max_error:
                    max_error = error
                    max_i1 = i1
            i1 = max_i1
            if i1 is not None:
                if self._take_step(i1, i2):
                    return 1
                
            i1 = random.randint(0, len(self.svm.alphas)-1)
            for i in range(len(self.svm.alphas)):
                if self._is_bound(self.svm.alphas[i1], self.C):
                    continue
                if self._take_step(i1, i2):
                    return 1
                i1 = (i1 + 1) % len(self.svm.alphas)

            i1 = random.randint(0, len(self.svm.alphas)-1)
            for i in range(len(self.svm.alphas)):
                if self._take_step(i1, i2):
                    return 1
                i1 = (i1 + 1) % len(self.svm.alphas)
        return 0

    def _take_step(self, i1, i2):
        if i1 == i2:
            return 0
        point = self.svm.xs
        alph1, alph2 = self.svm.alphas[i1], self.svm.alphas[i2]
        y1, y2 = self.svm.ys[i1], self.svm.ys[i2]
        E1 = self._E(i1)
        E2 = self._E(i2)
        s = y1*y2
        b = self.svm.b
        C = self.C
        eps = self.eps
        
        if y1 == y2:
            L = max(0.0, alph2+alph1-C)
            H = min(C, alph2+alph1)
        else:
            L = max(0.0, alph2-alph1)
            H = min(C, C+alph2-alph1)
        if L == H:
            return 0

        k11 = self.svm.kernel_fn(point[i1], point[i1])
        k12 = self.svm.kernel_fn(point[i1], point[i2])
        k22 = self.svm.kernel_fn(point[i2], point[i2])
        eta = k11+k22-2.0*k12
        if eta > 0.0:
            a2 = alph2 + y2*(E1-E2)/eta
            if a2 < L: a2 = L
            elif a2 > H: a2 = H
        else:
            f1 = y1*(E1+b) - alph1*k11 - s*alph2*k22
            f2 = y2*(E2+b) - s*alph1*k12 - alph2*k22
            L1 = alph1 + s*(alph2-L)
            H1 = alph1 + s*(alph2-H)
            Lobj = L1*f1 + L*f2 + 0.5*L1*L1*k11 + 0.5*L*L*k22 + s*L*L1*k12
            Hobj = H1*f1 + H*f2 + 0.5*H1*H1*k11 + 0.5*H*H*k22 + s*H*H1*k12
            
            if Lobj < Hobj-eps:
                a2 = L
            elif Lobj > Hobj+eps:
                a2 = H
            else:
                a2 = alph2

        if math.fabs(a2-alph2) < eps*(a2+alph2+eps):
            return 0

        a1 = alph1 + s*(alph2-a2)
        
        if not self._is_bound(a1, C):
            b1 = E1 + y1*(a1-alph1)*k11 + y2*(a2-alph2)*k12 + b
            b = b1
        elif not self._is_bound(a2, C):
            b2 = E2 + y1*(a1-alph1)*k12 + y2*(a2-alph2)*k22 + b
            b = b2
        else:
            b1 = E1 + y1*(a1-alph1)*k11 + y2*(a2-alph2)*k12 + b
            b2 = E2 + y1*(a1-alph1)*k12 + y2*(a2-alph2)*k22 + b
            b = (b1+b2)/2.0

        self.svm.alphas[i1], self.svm.alphas[i2] = a1, a2
        self.svm.b = b
        self._ecache = {}

        return 1

    def _E(self, i):
        if not self._ecache.has_key(i):
            self._ecache[i] = classify(self.svm, self.svm.xs[i])-self.svm.ys[i]
        return self._ecache[i]

    def _is_bound(self, alpha, C):
        return alpha <= 0 or alpha >= C

class KeerthiTrainer:
    """Sequential Minimal Optimization Trainer

    This is an implementation of Sequential Minimal Optimization plus
    the 2 modifications suggested by S.S. Keerthi, et al.:
    http://guppy.mpe.nus.edu.sg/~mpessk/

    Methods:
    train     Train a new SVM.

    """
    def train(self, training_set, results, kernel_fn, C, epsilon,
              update_fn=None):
        """train(self, training_set, results, kernel_fn, C, epsilon,
        update_fn=None) -> SVM

        """
        # Initializing variables
        self.svm = SVM(training_set, results,
                       [0.0]*len(training_set), 0.0, kernel_fn)
        self.eps = epsilon
        self.C = C
        self._fcache = {}

        # Initialize values.
        # b_low and b_up should be set to extreme values.
        # self.i_up and self.i_low should be set to a positive example,
        # and a negative example.
        self.b_low, self.b_up = 1, -1
        self.i_up = self.i_low = None
        i = 0
        while i < len(self.svm.alphas):
            if self.svm.ys[i] == 1 and self.i_up is None:
                self.i_up = i
            elif self.svm.ys[i] == -1 and self.i_low is None:
                self.i_low = i
            if self.i_up is not None and self.i_low is not None:
                break
            i = i + 1
        else:
            raise ValueError, \
                  "I could not find positive and negative training examples"

        # Training algorithm starts here.
        while 1:
            num_changed = 0
            for i2 in range(len(self.svm.alphas)):
                num_changed = num_changed + self._examine_example(i2)
            if not num_changed:
                break
            
            inner_loop_success = 1
            while (self.b_up <= self.b_low-2*self.eps) and \
                  inner_loop_success:
                i2 = self.i_low
                y2 = self.svm.ys[i2]
                alph2 = self.svm.alphas[i2]
                F2 = self._F(i2)
                i1 = self.i_up
                inner_loop_success = self._take_step(self.i_up, self.i_low)
                num_changed = num_changed + inner_loop_success
                
            if update_fn is not None:
                update_fn(num_changed, self.svm)

        # b_low and b_up define the range of allowable b's, within
        # a tolerance of epsilon.  I'll just choose the b in the middle.
        self.svm.b = (self.b_low+self.b_up)/2.0

        return self.svm

    def _examine_example(self, i2):
        y2 = self.svm.ys[i2]
        alph2 = self.svm.alphas[i2]

        F2 = self._F(i2)
        if self._in_I_0(i2):
            pass
        elif (self._in_I_1(i2) or self._in_I_2(i2)) and F2 < self.b_up:
            self.b_up = F2
            self.i_up = i2
        elif (self._in_I_3(i2) or self._in_I_4(i2)) and F2 > self.b_low:
            self.b_low = F2
            self.i_low = i2

        optimality = 1
        if self._in_I_0(i2) or self._in_I_1(i2) or self._in_I_2(i2):
            #if self.b_low-F2 > 2*self.eps:
            if self.b_low > F2:
                optimality = 0
                i1 = self.i_low
        if self._in_I_0(i2) or self._in_I_3(i2) or self._in_I_4(i2):
            #if F2-self.b_up > 2*self.eps:
            if F2 > self.b_up:
                optimality = 0
                i1 = self.i_up
        if optimality:
            return 0

        # for i2 in I_0 choose the better i1
        if self._in_I_0(i2):
            if self.b_low-F2 > F2-self.b_up:
                i1 = self.i_low
            else:
                i1 = self.i_up
        return self._take_step(i1, i2)

    def _take_step(self, i1, i2):
        if i1 == i2:
            return 0
        point = self.svm.xs
        alph1, alph2 = self.svm.alphas[i1], self.svm.alphas[i2]
        y1, y2 = self.svm.ys[i1], self.svm.ys[i2]
        F1, F2 = self._F(i1), self._F(i2)
        s = y1*y2
        b = self.svm.b
        C = self.C
        eps = self.eps
        
        if y1 == y2:
            L = max(0.0, alph2+alph1-C)
            H = min(C, alph2+alph1)
        else:
            L = max(0.0, alph2-alph1)
            H = min(C, C+alph2-alph1)
        if L == H:
            return 0

        k11 = self.svm.kernel_fn(point[i1], point[i1])
        k12 = self.svm.kernel_fn(point[i1], point[i2])
        k22 = self.svm.kernel_fn(point[i2], point[i2])
        eta = k11+k22-2.0*k12
        if eta > 0.0:
            a2 = alph2 + y2*(F1-F2)/eta
            if a2 < L: a2 = L
            elif a2 > H: a2 = H
        else:
            f1 = y1*F1 - alph1*k11 - s*alph2*k22
            f2 = y2*F2 - s*alph1*k12 - alph2*k22
            L1 = alph1 + s*(alph2-L)
            H1 = alph1 + s*(alph2-H)
            Lobj = L1*f1 + L*f2 + 0.5*L1*L1*k11 + 0.5*L*L*k22 + s*L*L1*k12
            Hobj = H1*f1 + H*f2 + 0.5*H1*H1*k11 + 0.5*H*H*k22 + s*H*H1*k12
            
            if Lobj < Hobj-eps:
                a2 = L
            elif Lobj > Hobj+eps:
                a2 = H
            else:
                a2 = alph2

        if math.fabs(a2-alph2) < eps*(a2+alph2+eps):
            return 0

        a1 = alph1 + s*(alph2-a2)
        
        self.svm.alphas[i1], self.svm.alphas[i2] = a1, a2
        self._fcache = {}

        i_low = i_up = b_low = b_up = None
        for i in range(len(self.svm.alphas)):
            if self._in_I_0(i) or self._in_I_3(i) or self._in_I_4(i):
                f = self._F(i)
                if b_low is None or f > b_low:
                    b_low = f
                    i_low = i
            if self._in_I_0(i) or self._in_I_1(i) or self._in_I_2(i):
                f = self._F(i)
                if b_up is None or f < b_up:
                    b_up = f
                    i_up = i
        self.i_low, self.i_up = i_low, i_up
        self.b_low, self.b_up = b_low, b_up

        return 1

    def _F(self, i):
        if not self._fcache.has_key(i):
            # assume b == 0
            self._fcache[i] = classify(self.svm, self.svm.xs[i])-self.svm.ys[i]
        return self._fcache[i]

    def _in_I_0(self, i):
        return self.svm.alphas[i] > 0 and self.svm.alphas[i] < self.C

    def _in_I_1(self, i):
        return self.svm.alphas[i] == 0 and self.svm.ys[i] == 1

    def _in_I_2(self, i):
        return self.svm.alphas[i] == self.C and self.svm.ys[i] == -1

    def _in_I_3(self, i):
        return self.svm.alphas[i] == self.C and self.svm.ys[i] == 1

    def _in_I_4(self, i):
        return self.svm.alphas[i] == 0 and self.svm.ys[i] == -1

class TransductiveTrainer:
    """Transductive Support Vector Machine Trainer

    This is used to train transductive support vector machines described
    in:
    Joachims, Thorsten.  Transductive Inference for Text Classification
    using Support Vector Machines.
    XXX url?

    The algorithm is modified from Platt's Sequential Minimal Optimization.

    Methods:
    train     Train a new SVM.

    """
    def train(self, training_set, results, num_test_cases,
              kernel_fn, C, C_negtest, C_postest, epsilon,
              update_fn=None):
        """train(self, training_set, results, num_test_cases,
        kernel_fn, C, C_negtest, C_postest, epsilon,
        update_fn=None) -> SVM

        """
        # Initializing variables
        self.svm = SVM(training_set, results,
                       [0.0]*len(training_set), 0.0, kernel_fn)
        self.eps = epsilon
        self.num_test_cases = num_test_cases
        self.C = C
        self.C_negtest, self.C_postest = C_negtest, C_postest
        self._ecache = {}

        # Training algorithm starts here.
        num_changed = 0
        examine_all = 1

        while num_changed > 0 or examine_all:
            num_changed = 0
            if examine_all:
                for i2 in range(len(self.svm.alphas)):
                    num_changed = num_changed + self._examine_example(i2)
            else:
                for i2 in range(len(self.svm.alphas)):
                    # skip the bound alphas
                    if self._is_bound(self.svm.alphas[i2], self._C(i2)):
                        continue
                    num_changed = num_changed + self._examine_example(i2)
            if update_fn is not None:
                update_fn(num_changed, self.svm)
            if examine_all:
                examine_all = 0
            elif num_changed == 0:
                examine_all = 1

        return self.svm

    def _examine_example(self, i2):
        y2 = self.svm.ys[i2]
        alph2 = self.svm.alphas[i2]
        E2 = self._E(i2)
        r2 = E2*y2

        if (r2 < -self.eps and alph2 < self._C(i2)) or \
           (r2 > self.eps and alph2 > 0):
            # sort it by E1 - E2
            max_i1 = None
            max_error = None
            for i1 in range(len(self.svm.alphas)):
                if i1 == i2 or self._is_bound(self.svm.alphas[i1],
                                              self._C(i1)):
                    continue
                E1 = self._E(i1)
                error = math.fabs(E1 - E2)
                if max_i1 is None or error > max_error:
                    max_error = error
                    max_i1 = i1
            i1 = max_i1
            if i1 is not None:
                if self._take_step(i1, i2):
                    return 1
                
            i1 = random.randint(0, len(self.svm.alphas)-1)
            for i in range(len(self.svm.alphas)):
                if self._is_bound(self.svm.alphas[i1], self._C(i1)):
                    continue
                if self._take_step(i1, i2):
                    return 1
                i1 = (i1 + 1) % len(self.svm.alphas)

            i1 = random.randint(0, len(self.svm.alphas)-1)
            for i in range(len(self.svm.alphas)):
                if self._take_step(i1, i2):
                    return 1
                i1 = (i1 + 1) % len(self.svm.alphas)
        return 0

    def _take_step(self, i1, i2):
        if i1 == i2:
            return 0
        point = self.svm.xs
        alph1, alph2 = self.svm.alphas[i1], self.svm.alphas[i2]
        y1, y2 = self.svm.ys[i1], self.svm.ys[i2]
        E1 = self._E(i1)
        E2 = self._E(i2)
        s = y1*y2
        b = self.svm.b
        eps = self.eps
        C1, C2 = self._C(i1), self._C(i2)

        if y1 == y2:
            L = max(0.0, alph2+alph1-C1)
            H = min(C2, alph2+alph1)
        else:
            L = max(0.0, alph2-alph1)
            H = min(C2, C1+alph2-alph1)
        if L == H:
            return 0

        k11 = self.svm.kernel_fn(point[i1], point[i1])
        k12 = self.svm.kernel_fn(point[i1], point[i2])
        k22 = self.svm.kernel_fn(point[i2], point[i2])
        eta = k11+k22-2.0*k12
        if eta > 0.0:
            a2 = alph2 + y2*(E1-E2)/eta
            if a2 < L: a2 = L
            elif a2 > H: a2 = H
        else:
            f1 = y1*(E1+b) - alph1*k11 - s*alph2*k22
            f2 = y2*(E2+b) - s*alph1*k12 - alph2*k22
            L1 = alph1 + s*(alph2-L)
            H1 = alph1 + s*(alph2-H)
            Lobj = L1*f1 + L*f2 + 0.5*L1*L1*k11 + 0.5*L*L*k22 + s*L*L1*k12
            Hobj = H1*f1 + H*f2 + 0.5*H1*H1*k11 + 0.5*H*H*k22 + s*H*H1*k12
            
            if Lobj < Hobj-eps:
                a2 = L
            elif Lobj > Hobj+eps:
                a2 = H
            else:
                a2 = alph2

        if math.fabs(a2-alph2) < eps*(a2+alph2+eps):
            return 0

        a1 = alph1 + s*(alph2-a2)

        if not self._is_bound(a1, C1):
            b1 = E1 + y1*(a1-alph1)*k11 + y2*(a2-alph2)*k12 + b
            b = b1
        elif not self._is_bound(a2, C2):
            b2 = E2 + y1*(a1-alph1)*k12 + y2*(a2-alph2)*k22 + b
            b = b2
        else:
            b1 = E1 + y1*(a1-alph1)*k11 + y2*(a2-alph2)*k12 + b
            b2 = E2 + y1*(a1-alph1)*k12 + y2*(a2-alph2)*k22 + b
            b = (b1+b2)/2.0

        self.svm.alphas[i1], self.svm.alphas[i2] = a1, a2
        self.svm.b = b
        self._ecache = {}

        return 1

    def _E(self, i):
        if not self._ecache.has_key(i):
            self._ecache[i] = classify(self.svm, self.svm.xs[i])-self.svm.ys[i]
        return self._ecache[i]

    def _C(self, i):
        if i < len(self.svm.ys)-self.num_test_cases:
            return self.C
        elif self.svm.ys[i] > 0:
            return self.C_postest
        return self.C_negtest

    def _is_bound(self, alpha, C):
        return alpha <= 0 or alpha >= C

def train(training_set, results,
          kernel_fn=LinearKernel(), C=1.0, epsilon=1E-3, update_fn=None):
    """train(training_set, results,
    kernel_fn=LinearKernel(), C=1.0, epsilon=1E-3, update_fn=None) -> SVM

    Train a new support vector machine.  training_set is a list
    containing a list of numbers, each which is a vector in the
    training set.  results is a parallel list that contains either
    1 or -1, which describes the class that the training vector belongs
    to.

    kernel_fn, C, and epsilon are optional parameters used to tune the
    support vector machine.

    update_fn is an optional callback function that called at the
    end of every round of training.  It should take 2 parameters:
        update_fn(num_changed, svm)

    """
    svm = KeerthiTrainer().train(
        training_set, results, kernel_fn, C, epsilon, update_fn=update_fn)
    return svm

def train_transductive(training_set, results, test_set, num_pos,
                       kernel_fn=LinearKernel(), C=1.0, C_test=0.5,
                       epsilon=1E-3, update_fn=None):
    """train_transductive(training_set, results, test_set, num_pos,
    kernel_fn=LinearKernel(), C=1.0, C_test=0.5,
    epsilon=1E-3, update_fn=None) -> SVM

    XXX document this

    """
    svm = train(training_set, results,
                kernel_fn=kernel_fn, C=C, epsilon=epsilon, update_fn=update_fn)

    # Assign the top num_pos results into the positive class.
    # Everything else is in the negative.
    test_scores = []   # list of tuples (score, index)
    for i in range(len(test_set)):
        score = classify(svm, test_set[i])
        test_scores.append((score, i))
    test_scores.sort()
    test_results = [None] * len(test_scores)
    for i in range(len(test_scores)):
        score, index = test_scores[i]
        if i < num_pos:
            test_results[index] = 1
        else:
            test_results[index] = -1
    
    trainer = TransductiveTrainer()
    C_neg = 1E-5
    C_pos = 1E-5 * num_pos/(len(test_set)-num_pos)
    while C_neg < C_test or C_pos < C_test:
        svm = trainer.train(training_set+test_set, results+test_results,
                            len(test_set), kernel_fn, C, C_neg, C_pos,
                            epsilon, update_fn=update_fn)
        changed = 1
        while changed:
            changed = 0
            # Look for an m and l that should be in different classes, and
            # that are both misclassified.  If I see such a case, then
            # switch them.
            #
            # The amount misclassified is:
            # gamma = 1 - y*score
            # If the sample is misclassified, y*score is negative.
            # If it's within the margin, y*score is positive and < 1.
            # Thus, I need to check to make sure gamma > 0, and the
            # sum of the gammas > 2.

            total = len(test_set)**2  # total number of combinations to try
            offset = random.randint(0, total-1)  # start at a random point
            for i in range(total):
                index = (i + offset) % total
                m = index / len(test_set)
                l = index % len(test_set)

                # Make sure m and l should be in different classes.
                if test_results[m] == test_results[l]:
                    continue
                # Make sure they're both misclassified.
                gamma_m = 1.0 - test_results[m]*classify(svm, test_set[m])
                gamma_l = 1.0 - test_results[l]*classify(svm, test_set[l])
                if gamma_m <= 0 or gamma_l <= 0 or (gamma_m+gamma_l) < 2:
                    continue
                
                # Switch them, and reclassify.
                test_results[l], test_results[m] = \
                                 test_results[m], test_results[l]
                
                svm = trainer.train(training_set+test_set,results+test_results,
                                    len(test_set), kernel_fn, C, C_neg, C_pos,
                                    epsilon, update_fn=update_fn)
                changed = 1
                break
                    
        C_neg = min(C_neg*2, C_test)
        C_pos = max(C_pos*2, C_test)
        
    return svm
    
def _dot(x, y):
    """_dot(x, y) -> x*y

    Return the dot product of x and y, where x and y are lists of numbers.

    """
    if len(x) != len(y):
        raise ValueError, "vectors must be same length"
    sum = 0.0
    for i in range(len(x)):
        sum = sum + x[i]*y[i]
    return sum

def _sparse_dot(x, y):
    sum = 0.0
    for k in x.keys():
        if y.has_key(k):
            sum = sum + x[k] * y[k]
    return sum

def _subtract(x, y):
    """_subtract(x, y) -> x-y

    Return x subtract y, where x and y are lists of numbers.

    """
    if len(x) != len(y):
        raise ValueError, "vectors must be same length"
    s = [None] * len(x)
    for i in range(len(x)):
        s[i] = x[i]-y[i]
    return s

def _sparse_subtract(x, y):
    done = {}  # have I handled this index yet?
    vec = {}
    for k in x.keys():
        if y.has_key(k):
            vec[k] = x[k] - y[k]
        else:
            vec[k] = x[k]
        done[k] = 1
    for k in y.keys():
        if done.has_key(k):
            continue
        vec[k] = -y[k]
    return vec

def _sign(x):
    """_sign(x) -> 1 or -1

    Return 1/-1 depending on the sign of x.

    """
    if x >= 0:
        return 1
    return -1

import cSVM
_sparse_dot = cSVM._sparse_dot
_dot = cSVM._dot
classify = cSVM.classify
