From 85a4962556b67d1cc0668ecb2fbb03b3b4dd6e7e Mon Sep 17 00:00:00 2001
From: Joe Zhao <ztuowen@gmail.com>
Date: Sat, 11 Apr 2015 17:46:19 +0800
Subject: completed & tested, train & predict

---
 model/ranksvmtn.cpp | 117 ++++++++++++++++++++++++++--------------------------
 1 file changed, 58 insertions(+), 59 deletions(-)

(limited to 'model/ranksvmtn.cpp')

diff --git a/model/ranksvmtn.cpp b/model/ranksvmtn.cpp
index 776d4db..959ea7d 100644
--- a/model/ranksvmtn.cpp
+++ b/model/ranksvmtn.cpp
@@ -7,7 +7,8 @@ using namespace std;
 using namespace Eigen;
 
 const int maxiter = 10;
-const double prec=1e-3;
+const double prec=1e-4;
+const double C=1;
 
 int cg_solve(const MatrixXd &A, const VectorXd &b, VectorXd &x)
 {
@@ -20,9 +21,7 @@ int cg_solve(const MatrixXd &A, const VectorXd &b, VectorXd &x)
     {
         // Non preconditioned version
         r_1 = res.dot(res);
-        cout<<step<<":"<<r_1<<endl;
-        write_stream(cout,res);
-        if (r_1<1e-5) // Terminate condition
+        if (r_1<1e-10) // Terminate condition
             break;
         if (step){
             beta = r_1 / r_2;
@@ -33,9 +32,6 @@ int cg_solve(const MatrixXd &A, const VectorXd &b, VectorXd &x)
         alpha = r_1/p.dot(q);
         x=x+p*alpha;
         res=res-q*alpha;
-        write_stream(cout,p);
-        write_stream(cout,q);
-        cin.get();
         ++step;
         r_2=r_1;
     }
@@ -43,16 +39,12 @@ int cg_solve(const MatrixXd &A, const VectorXd &b, VectorXd &x)
 }
 
 // Calculate objfunc gradient & support vectors
-int objfunc_linear(const VectorXd &w,const MatrixXd &A,const double C,VectorXd &pred,VectorXd &grad, double &obj,MatrixXd &sv)
+int objfunc_linear(const VectorXd &w,const MatrixXd &D,const MatrixXd &A,const double C,VectorXd &pred,VectorXd &grad, double &obj)
 {
-    pred = pred.cwiseMax(MatrixXd::Zero(pred.rows(),pred.cols()));
+    for (int i=0;i<pred.rows();++i)
+        pred(i)=pred(i)>0?pred(i):0;
     obj = (pred.cwiseProduct(pred)*C).sum()/2 + w.dot(w)/2;
-    grad = w - (((pred*C).transpose()*A)*w).transpose();
-    for (int i=0;i<pred.cols();++i)
-        if (pred(i)>0)
-            sv(i,i)=1;
-        else
-            sv(i,i)=0;
+    grad = w - (((pred*C).transpose()*A)*D).transpose();
     return 0;
 }
 
@@ -63,36 +55,40 @@ int line_search(const VectorXd &w,const MatrixXd &D,const MatrixXd &A,const Vect
     double g,h;
     t = 0;
     VectorXd Xd=A*(D*step);
+    VectorXd pred2;
     while (1)
     {
-        pred = pred - t*Xd;
+        pred2 = pred - t*Xd;
         g=wd+t*dd;
         h=dd;
-        for (int i=0;i<pred.cols();++i)
-            if (pred(i)>0) {
-                g += pred(i)*Xd(i);
-                h += Xd(i)*Xd(i);
+        for (int i=0;i<pred2.rows();++i)
+            if (pred2(i)>0) {
+                g -= C*pred2(i)*Xd(i);
+                h += C*Xd(i)*Xd(i);
             }
+        g=g+1e-12;
+        h=h+1e-12;
+        t=t-g/h;
+        cout<<g<<":"<<h<<endl;
+        cin.get();
         if (g*g/h<1e-10)
             break;
     }
+    pred=pred2;
     return 0;
 }
 
 int train_orig(int fsize, MatrixXd &D,MatrixXd &A,VectorXd &weight){
     int iter = 0;
-    double C=1;
 
-    long n=D.rows();
-    LOG(INFO) << "training with feature size:" << fsize << " Data size:" << n;
-    MatrixXd sv=MatrixXd::Identity(n, n);
+    long n=A.rows();
+    LOG(INFO) << "training with feature size:" << fsize << " Data size:" << n << " Relation size:" << A.rows();
     VectorXd grad(fsize);
     VectorXd step(fsize);
     VectorXd pred(n);
     double obj,t;
 
     pred=VectorXd::Ones(n) - (A*(D*weight));
-
     while (true)
     {
         iter+=1;
@@ -103,21 +99,28 @@ int train_orig(int fsize, MatrixXd &D,MatrixXd &A,VectorXd &weight){
         }
 
         // Generate support vector matrix sv & gradient
-        objfunc_linear(weight,A,C,pred,grad,obj,sv);
+        objfunc_linear(weight,D,A,C,pred,grad,obj);
         step = grad*0;
-        MatrixXd H = MatrixXd::Identity(grad.cols(),grad.cols());
-
+        MatrixXd H = MatrixXd::Identity(grad.rows(),grad.rows());
         // Compute Hessian directly
         for (int i=0;i<n;++i)
-            if (sv(i,i)>0)
-            H = H + 2*C*A.row(i).transpose()*A.row(i);
+            if (pred(i)>0) {
+                VectorXd v = A.row(i)*D;
+                H = H + C * (v * v.transpose());
+            }
         // Solve
+        //cout<<obj<<endl;
         cg_solve(H,grad,step);
         // do line search
         line_search(weight,D,A,step,pred,t);
         weight=weight+step*t;
+        int sv=0;
+        for (int i=0;i<n;++i)
+            if (pred(i)>0)
+                ++sv;
         // When dec is small enough
-        if (-step.dot(grad) < prec * obj)
+        LOG(INFO)<<"Iter: "<<iter<<" Obj: " <<obj<<" SV: "<< sv << " Newton decr:"<<step.dot(grad)/2 << " linesearch: "<< -t ;
+        if (step.dot(grad) < prec * obj)
             break;
     }
     return 0;
@@ -125,55 +128,51 @@ int train_orig(int fsize, MatrixXd &D,MatrixXd &A,VectorXd &weight){
 
 int RSVMTN::train(DataList &D){
     MatrixXd Data(D.getSize(),D.getfSize()),A;
-    int i=0,j=0;
-    list<DataEntry*>::iterator iter,st,nx;
-    for (iter= D.getData().begin();i<D.getSize();++i,++iter)
-        for (j=0;j<D.getfSize();++j)
-            Data(i,j)=(*iter)->feature(j);
-    nx=st=iter= D.getData().begin();
-    ++nx;
+    int i,j;
+    LOG(INFO)<<"Processing input";
+    for (i=0;i<D.getSize();++i) {
+        for (j = 0; j < D.getfSize(); ++j)
+            Data(i, j) = (D.getData()[i])->feature(j);
+    }
     int cnt=0;
-    while (iter!=D.getData().end())
+    i=j=0;
+    while (i<D.getSize())
     {
-        if ((nx == D.getData().end())||(*iter)->qid!=(*nx)->qid)
+        if ((i+1 == D.getSize())|| D.getData()[i]->qid!=D.getData()[i+1]->qid)
         {
-            list<DataEntry*>::iterator high,low=iter;
-            for (high=st;((*high)->rank)>0;++high)
-                for (low=iter;((*low)->rank)<0;--low)
-                    ++cnt;
-            st = nx;
+            int high=j;
+            while (D.getData()[high]->rank>0)
+                ++high;
+            cnt += (high-j)*(i-high+1);
+            j = i+1;
         }
-        ++iter;
+        ++i;
     }
     A.resize(cnt,D.getSize());
-    nx=st=iter= D.getData().begin();
-    ++nx;
     cnt=i=j=0;
-    while (iter!=D.getData().end())
+    while (i<D.getSize())
     {
-        if ((nx == D.getData().end())||(*iter)->qid!=(*nx)->qid)
+        if ((i+1 == D.getSize())|| D.getData()[i]->qid!=D.getData()[i+1]->qid)
         {
             int v1=j,v2;
-            list<DataEntry*>::iterator high,low=iter;
-            for (high=st;((*high)->rank)>0;++high,++v1)
-                for (low=iter,v2=i;((*low)->rank)<0;--low,--v2) {
+            for (v1=j;(D.getData()[v1]->rank)>0;++v1)
+                for (v2=i;(D.getData()[v2]->rank)<0;--v2) {
                     A(cnt,v1) = 1;
                     A(cnt,v2) = -1;
                     ++cnt;
                 }
-            st = nx;
-            j=i+1;
+            j = i+1;
         }
         ++i;
-        ++iter;
     }
     train_orig(fsize,Data,A,model.weight);
     return 0;
 };
 
-int RSVMTN::predict(DataList &D, list<double> &res){
+int RSVMTN::predict(DataList &D, vector<double> &res){
     //TODO define A
-    for (list<DataEntry*>::iterator i=D.getData().begin(), end=D.getData().end();i!=end;++i)
-        res.push_back(((*i)->feature).dot(model.weight));
+    res.clear();
+    for (int i=0;i<D.getSize();++i)
+        res.push_back(((D.getData()[i])->feature).dot(model.weight));
     return 0;
 };
\ No newline at end of file
-- 
cgit v1.2.3-70-g09d2