inf.small<-10^(-6)
K<-function(x,y,sigma2=ifelse(is.vector(x),1,ncol(x)),kernel="radial") {
   if(is.vector(x)) {
       if (kernel=="linear") { return ( sum(x*y) ) } else { return ( exp(-sum((x-y)^2)/sigma2) ) }
   } else {
       if (kernel=="linear") return ( x %*% t(y) ) else {
          res<-matrix(0,nrow(x),nrow(y))
          for ( i in 1:nrow(x) ) {
          for ( j in 1:nrow(y) ) { 
             res[i,j]<-exp(-sum((x[i,]-y[j,])^2)/sigma2) 
          }
          }
          return (res)
       }
   }
}

tsvm.dca<-function(data.l,data.u,cost1=1,cost2=1,kernel="radial") {
   require(quadprog)	
   require(lpSolve)	
   require(e1071)
	
   data.l<-data.matrix(data.l)
   x.l<-data.l[,-1]
   y.l<-data.l[,1]
   x.u<-data.matrix(data.u)
   p<-ncol(x.l)
   l<-nrow(x.l)
   u<-nrow(x.u)

   m.ini<-svm(x.l,y.l,cost=cost1,kernel="radial")		#Initialization
   alpha.ini<-rep(0,length(y.l))
   alpha.ini[m.ini$index]<-m.ini$coefs
   alpha.ini<-alpha.ini*y.l
   coef<-as.vector(apply(as.vector(alpha.ini)*y.l*x.l,2,sum))
   const<--m.ini$rho
   grad<-cost2/u*ifelse(as.vector(x.u%*%coef)+const>0,1,-1)
   fun<-cost1/l*sum(pos(1-y.l*(as.vector(x.l%*%coef)+const)))+cost2/u*sum(pos(1-abs(as.vector(x.u%*%coef)+const)))+sum(coef^2)/2
   f.u<-as.vector(x.u%*%coef)
   f.l<-as.vector(x.l%*%coef)

   fun.old<-0
   fun.vec<-0.01
   epsilon<-10^-2
   while(abs(fun-fun.old)>epsilon*fun & abs(sum(grad))<cost2 & sum(abs(fun-fun.vec)<(epsilon*fun))<1  ) {
      coef.old<-coef
      const.old<-const
      fun.old<-fun
#      print(fun)
      fun.vec<-c(fun.vec,fun)

      Dt<-cbind(rbind(diag(rep(1,l)),matrix(0,u,l)),rbind(matrix(0,l,u),diag(rep(1,u))),rbind(matrix(0,l,u),diag(rep(-1,u))))
      Dmat1<-K(y.l*x.l,y.l*x.l)
      Dmat2<--1*K(y.l*x.l,x.u) 
      Dmat3<-K(x.u,x.u) 
      Dmat<-t(Dt) %*% cbind(rbind(Dmat1,t(Dmat2)),rbind(Dmat2,Dmat3)) %*% Dt  + diag(rep(inf.small,l+2*u))
      dvec<--1*c(1-y.l*(K(x.l,x.u)%*%grad), -1+K(x.u,x.u)%*%grad, -1-K(x.u,x.u)%*%grad)
      Amat<-t(rbind(t(c(y.l,rep(-1,u)))%*%Dt, cbind(diag(rep(1,l)),matrix(0,l,2*u)), cbind(diag(rep(-1,l)),matrix(0,l,2*u)), cbind(matrix(0,u,l),diag(rep(1,u)),matrix(0,u,u)), cbind(matrix(0,u,l+u),diag(rep(1,u))), cbind(matrix(0,u,l),diag(rep(-1,u)),diag(rep(-1,u)))))
      bvec<-c(-sum(grad),rep(0,l),rep(-cost1/l,l),rep(0,u),rep(0,u),rep(-cost2/u,u))
      meq<-1

      sol<-solve.QP(Dmat, dvec, Amat, bvec, meq)$solution
      if (as.character(sum(sol))=="NaN") break
      f.l<-as.vector(apply(sol[1:l]*y.l*K(x.l,x.l),2,sum) + apply(as.vector(cbind(diag(rep(1,u)), diag(rep(-1,u))) %*% sol[l+(1:(2*u))])*K(x.u,x.l),2,sum) + apply(grad*K(x.u,x.l),2,sum))
      f.u<-as.vector(apply(sol[1:l]*y.l*K(x.l,x.u),2,sum) + apply(as.vector(cbind(diag(rep(1,u)), diag(rep(-1,u))) %*% sol[l+(1:(2*u))])*K(x.u,x.u),2,sum) + apply(grad*K(x.u,x.u),2,sum))

      i0<-order(abs(sol[1:l]-cost1/2/l))[1]
      if ( abs(sol[i0]-cost1/l/2)-cost1/l/2>10^-6) {
         const<-1/y.l[i0] - f.l[i0,]
      } else {
         f.obj<-c(-sum(grad),rep(cost1/l,l),rep(cost2/u,u))
         f.con<-rbind(cbind(y.l,diag(rep(1,l)),matrix(0,l,u)), cbind(rep(0,l),diag(rep(1,l)), matrix(0,l,u)), cbind(rep(-1,u), matrix(0,u,l),diag(rep(1,u))), cbind(rep(1,u),matrix(0,u,l),diag(rep(1,u))), cbind(matrix(0,u,l+1),diag(1,u,u)))
         f.dir<-rep(">=",2*l+3*u)
         f.rhs<-c(1-y.l*f.l,rep(0,l),-1+f.u,-1-f.u,rep(0,u))
         const.sol<-lp("min", f.obj, f.con, f.dir, f.rhs)$solution
         const<-const.sol[1]
      }

      grad<-cost2/u*ifelse(as.vector(f.u+const)>0,1,-1)

      fun<-as.numeric(cost1/l*sum(pos(1-y.l*as.vector(f.l+const)))+cost2/u*sum(pos(1-abs(as.vector(f.u+const))))+0.5 * t(sol) %*% Dmat %*% sol)
   }

   return(list(coef.l=sol[1:l]*y.l, coef.u=as.vector(cbind(diag(rep(1,u)),diag(rep(-1,u))) %*% sol[l+(1:(2*u))]), gradient=grad, const=const, class.u=sign(f.u+const)))
}