# _*_ coding: utf-8 _*_
from numpy import *

## SVM适用于二分类（数值型和标称型数据）
## 道理很深奥
def kernelTrans(X, A, kTup): #calc the kernel or transform data to a higher dimensional space
	m, n = shape(X)
	K = mat(zeros((m, 1)))
	if kTup[0] == 'lin':
		K = X * A.T   #linear kernel
	elif kTup[0] == 'rbf':
		for j in range(m):
			deltaRow = X[j, :] - A
			K[j] = deltaRow * deltaRow.T
		K = exp(K / (-1 * kTup[1] ** 2)) #divide in NumPy is element-wise not matrix like Matlab
	else:
		raise NameError('Houston We Have a Problem -- That Kernel is not recognized')
	return K


class RBF:
	def __init__(self, svInd, sVs, labelSV, alphas, b):
		self.svInd = svInd
		self.sVs = sVs
		self.labelSV = labelSV
		self.alphas = alphas
		self.b = b


class SVMLib:
	def selectJrand(self, i, m):
		j = i #we want to select any J not equal to i
		while (j == i):
			j = int(random.uniform(0, m))
		return j

	def clipAlpha(self, aj, H, L):
		if aj > H:
			aj = H
		if L > aj:
			aj = L
		return aj

	def smoSimple(self, dataMatIn, classLabels, C, toler, maxIter):
		dataMatrix = mat(dataMatIn);
		labelMat = mat(classLabels).transpose()
		b = 0;
		m, n = shape(dataMatrix)
		alphas = mat(zeros((m, 1)))
		iter = 0
		while (iter < maxIter):
			alphaPairsChanged = 0
			for i in range(m):
				fXi = float(multiply(alphas, labelMat).T * (dataMatrix * dataMatrix[i, :].T)) + b
				Ei = fXi - float(labelMat[i])#if checks if an example violates KKT conditions
				if ((labelMat[i] * Ei < -toler) and (alphas[i] < C)) or (
					(labelMat[i] * Ei > toler) and (alphas[i] > 0)):
					j = self.selectJrand(i, m)
					fXj = float(multiply(alphas, labelMat).T * (dataMatrix * dataMatrix[j, :].T)) + b
					Ej = fXj - float(labelMat[j])
					alphaIold = alphas[i].copy();
					alphaJold = alphas[j].copy();
					if (labelMat[i] != labelMat[j]):
						L = max(0, alphas[j] - alphas[i])
						H = min(C, C + alphas[j] - alphas[i])
					else:
						L = max(0, alphas[j] + alphas[i] - C)
						H = min(C, alphas[j] + alphas[i])
					if L == H: print "L==H"; continue
					eta = 2.0 * dataMatrix[i, :] * dataMatrix[j, :].T - dataMatrix[i, :] * dataMatrix[i,
					                                                                       :].T - dataMatrix[j,
					                                                                              :] * dataMatrix[j,
					                                                                                   :].T
					if eta >= 0: print "eta>=0"; continue
					alphas[j] -= labelMat[j] * (Ei - Ej) / eta
					alphas[j] = self.clipAlpha(alphas[j], H, L)
					if (abs(alphas[j] - alphaJold) < 0.00001): print "j not moving enough"; continue
					alphas[i] += labelMat[j] * labelMat[i] * (alphaJold - alphas[j])#update i by the same amount as j
					#the update is in the oppostie direction
					b1 = b - Ei - labelMat[i] * (alphas[i] - alphaIold) * dataMatrix[i, :] * dataMatrix[i, :].T - \
					     labelMat[j] * (alphas[j] - alphaJold) * dataMatrix[i, :] * dataMatrix[j, :].T
					b2 = b - Ej - labelMat[i] * (alphas[i] - alphaIold) * dataMatrix[i, :] * dataMatrix[j, :].T - \
					     labelMat[j] * (alphas[j] - alphaJold) * dataMatrix[j, :] * dataMatrix[j, :].T
					if (0 < alphas[i]) and (C > alphas[i]):
						b = b1
					elif (0 < alphas[j]) and (C > alphas[j]):
						b = b2
					else:
						b = (b1 + b2) / 2.0
					alphaPairsChanged += 1
					print "iter: %d i:%d, pairs changed %d" % (iter, i, alphaPairsChanged)
			if (alphaPairsChanged == 0):
				iter += 1
			else:
				iter = 0
			print "iteration number: %d" % iter
		return b, alphas

	def kernelTrans(X, A, kTup): #calc the kernel or transform data to a higher dimensional space
		m, n = shape(X)
		K = mat(zeros((m, 1)))
		if kTup[0] == 'lin':
			K = X * A.T   #linear kernel
		elif kTup[0] == 'rbf':
			for j in range(m):
				deltaRow = X[j, :] - A
				K[j] = deltaRow * deltaRow.T
			K = exp(K / (-1 * kTup[1] ** 2)) #divide in NumPy is element-wise not matrix like Matlab
		else:
			raise NameError('Houston We Have a Problem -- \
		That Kernel is not recognized')
		return K

	class optStruct:
		def __init__(self, dataMatIn, classLabels, C, toler, kTup):  # Initialize the structure with the parameters
			self.X = dataMatIn
			self.labelMat = classLabels
			self.C = C
			self.tol = toler
			self.m = shape(dataMatIn)[0]
			self.alphas = mat(zeros((self.m, 1)))
			self.b = 0
			self.eCache = mat(zeros((self.m, 2))) #first column is valid flag
			self.K = mat(zeros((self.m, self.m)))
			for i in range(self.m):
				self.K[:, i] = kernelTrans(self.X, self.X[i, :], kTup)

	def calcEk(self, oS, k):
		fXk = float(multiply(oS.alphas, oS.labelMat).T * oS.K[:, k] + oS.b)
		Ek = fXk - float(oS.labelMat[k])
		return Ek

	def selectJ(self, i, oS, Ei):         #this is the second choice -heurstic, and calcs Ej
		maxK = -1;
		maxDeltaE = 0;
		Ej = 0
		oS.eCache[i] = [1, Ei]  #set valid #choose the alpha that gives the maximum delta E
		validEcacheList = nonzero(oS.eCache[:, 0].A)[0]
		if (len(validEcacheList)) > 1:
			for k in validEcacheList:   #loop through valid Ecache values and find the one that maximizes delta E
				if k == i: continue #don't calc for i, waste of time
				Ek = self.calcEk(oS, k)
				deltaE = abs(Ei - Ek)
				if (deltaE > maxDeltaE):
					maxK = k;
					maxDeltaE = deltaE;
					Ej = Ek
			return maxK, Ej
		else:   #in this case (first time around) we don't have any valid eCache values
			j = self.selectJrand(i, oS.m)
			Ej = self.calcEk(oS, j)
		return j, Ej

	def updateEk(self, oS, k):#after any alpha has changed update the new value in the cache
		Ek = self.calcEk(oS, k)
		oS.eCache[k] = [1, Ek]

	def innerL(self, i, oS):
		Ei = self.calcEk(oS, i)
		if ((oS.labelMat[i] * Ei < -oS.tol) and (oS.alphas[i] < oS.C)) or (
			(oS.labelMat[i] * Ei > oS.tol) and (oS.alphas[i] > 0)):
			j, Ej = self.selectJ(i, oS, Ei) #this has been changed from selectJrand
			alphaIold = oS.alphas[i].copy();
			alphaJold = oS.alphas[j].copy();
			if (oS.labelMat[i] != oS.labelMat[j]):
				L = max(0, oS.alphas[j] - oS.alphas[i])
				H = min(oS.C, oS.C + oS.alphas[j] - oS.alphas[i])
			else:
				L = max(0, oS.alphas[j] + oS.alphas[i] - oS.C)
				H = min(oS.C, oS.alphas[j] + oS.alphas[i])
			if L == H: print "L==H"; return 0
			eta = 2.0 * oS.K[i, j] - oS.K[i, i] - oS.K[j, j] #changed for kernel
			if eta >= 0: print "eta>=0"; return 0
			oS.alphas[j] -= oS.labelMat[j] * (Ei - Ej) / eta
			oS.alphas[j] = self.clipAlpha(oS.alphas[j], H, L)
			self.updateEk(oS, j) #added this for the Ecache
			if (abs(oS.alphas[j] - alphaJold) < 0.00001): print "j not moving enough"; return 0
			oS.alphas[i] += oS.labelMat[j] * oS.labelMat[i] * (
			alphaJold - oS.alphas[j])#update i by the same amount as j
			self.updateEk(oS, i) #added this for the Ecache                    #the update is in the oppostie direction
			b1 = oS.b - Ei - oS.labelMat[i] * (oS.alphas[i] - alphaIold) * oS.K[i, i] - oS.labelMat[j] * (
			oS.alphas[j] - alphaJold) * oS.K[i, j]
			b2 = oS.b - Ej - oS.labelMat[i] * (oS.alphas[i] - alphaIold) * oS.K[i, j] - oS.labelMat[j] * (
			oS.alphas[j] - alphaJold) * oS.K[j, j]
			if (0 < oS.alphas[i]) and (oS.C > oS.alphas[i]):
				oS.b = b1
			elif (0 < oS.alphas[j]) and (oS.C > oS.alphas[j]):
				oS.b = b2
			else:
				oS.b = (b1 + b2) / 2.0
			return 1
		else:
			return 0

	def smoP(self, dataMatIn, classLabels, C, toler, maxIter, kTup=('lin', 0)):    #full Platt SMO
		oS = self.optStruct(mat(dataMatIn), mat(classLabels).transpose(), C, toler, kTup)
		iter = 0
		entireSet = True;
		alphaPairsChanged = 0
		while (iter < maxIter) and ((alphaPairsChanged > 0) or (entireSet)):
			alphaPairsChanged = 0
			if entireSet:   #go over all
				for i in range(oS.m):
					alphaPairsChanged += self.innerL(i, oS)
					print "fullSet, iter: %d i:%d, pairs changed %d" % (iter, i, alphaPairsChanged)
				iter += 1
			else:#go over non-bound (railed) alphas
				nonBoundIs = nonzero((oS.alphas.A > 0) * (oS.alphas.A < C))[0]
				for i in nonBoundIs:
					alphaPairsChanged += self.innerL(i, oS)
					print "non-bound, iter: %d i:%d, pairs changed %d" % (iter, i, alphaPairsChanged)
				iter += 1
			if entireSet:
				entireSet = False #toggle entire set loop
			elif (alphaPairsChanged == 0):
				entireSet = True
			print "iteration number: %d" % iter
		return oS.b, oS.alphas

	def calcWs(self, alphas, dataArr, classLabels):
		X = mat(dataArr);
		labelMat = mat(classLabels).transpose()
		m, n = shape(X)
		w = zeros((n, 1))
		for i in range(m):
			w += multiply(alphas[i] * labelMat[i], X[i, :].T)
		return w

	def fit(self, dataMatIn, classLabels, C, toler, maxIter, kTup=('lin', 0)):
		b, alpha = self.smoP(dataMatIn, classLabels, C, toler, maxIter, kTup=('lin', 0))
		return self.calcWs(alpha, dataMatIn, classLabels), b

	def predict(self, data, w, b):
		return data * mat(w) + b

	def fit_RBF(self, dataArr, labelArr, C, toler, maxIter, kTup):
		b, alphas = self.smoP(dataArr, labelArr, C, toler, maxIter, kTup)
		#C=200 important
		dataMatIn = mat(dataArr);
		svInd = nonzero(alphas.A > 0)[0]
		sVs = dataMatIn[svInd] #get matrix of only support vectors
		labelMat = mat(labelArr).transpose()
		labelSV = labelMat[svInd];

		return RBF(svInd, sVs, labelSV, alphas, b)

	def predict_RBF(self, RBF, datMat, kTup):
		kernelEval = kernelTrans(RBF.sVs, datMat, kTup)
		predict = kernelEval.T * multiply(RBF.labelSV, RBF.alphas[RBF.svInd]) + RBF.b
		return predict

