PLearn 0.1
linearalign.h
Go to the documentation of this file.
00001 
00002 #ifndef linearalign_INC
00003 #define linearalign_INC
00004 
00005 
00006 #include <plearn/io/load_and_save.h>
00007 #include <plearn/math/TMat_maths.h>
00008 #include "SurfaceMesh.h"
00009 #include "geometry.h"
00010 #include <plearn/misc/qld_interface.h>
00011 #include <plearn/ker/GaussianKernel.h>
00012 //#include "Molecule.h"
00013 #include "Template.h"
00014 extern "C"{
00015 #include <glpk.h>
00016 }
00017 using namespace PLearn;
00018 using namespace std;
00019 
00026 void static nodekernel(const Mat& x,const Mat& y,const Mat& dev,real sigma,Mat& ans){
00027         int xrows=x.nrows();
00028         int yrows=y.nrows();
00029         int cols=x.ncols();
00030         if(cols!=y.ncols()){
00031                 PLERROR("The property matrices have different number of fields");
00032         }
00033         ans=Mat(xrows,yrows,0.0);
00034         //cout<<"number of columns "<<cols<<endl;
00035         for(int i=0;i<xrows;i++){
00036                 for(int j=0;j<yrows;j++){
00037                         for(int k=0;k<x.ncols();k++){
00038                                 real term = (x[i][k]-y[j][k])/(dev[j][k]*sigma);
00039                                 term *= term;
00040                                 ans[i][j] += term;
00041                         }
00042                         ans[i][j] = exp(-ans[i][j]);
00043                 }
00044         }
00045 }
00046 
00047 bool static compare(const pair<real,int>& a,const pair<real,int>& b){
00048         if (a.first < b.first) return true;
00049         else return false;
00050 }
00051 
00055 void static sortedIndexList(const Vec& tosort,vector<int>& slist){
00056         int n = tosort.length();
00057         vector< pair<real,int> > ilist;
00058         for(int i=0;i<n;i++){
00059                 pair<real,int> temp;
00060                 temp.first = tosort[i];
00061                 temp.second = i;
00062                 ilist.push_back(temp);
00063         }
00064         sort(ilist.begin(),ilist.end(),compare);
00065         slist = vector<int>(n);
00066         for(int i=0;i<n;i++){
00067                 slist[i] = ilist[i].second;
00068         }
00069 }
00070 
00077 void static findRelevantWeights(const Mat& nkmat,vector< vector<bool> >& weights , int n){
00078         const int nx = nkmat.nrows();
00079         const int ny = nkmat.ncols();
00080         weights = vector< vector <bool> >(nx);
00081         for(int i=0;i<nx;i++){
00082                 weights[i] = vector<bool>(ny,false);
00083         }       
00084         /*n weights in each row be present*/
00085         
00086         for(int i=0;i<nx;i++){
00087                 vector< pair<real,int> > w(ny);
00088                 for(int j=0;j<ny;j++){
00089                         w[j].first=nkmat[i][j];
00090                         w[j].second=j;
00091                 }
00092                 sort(w.begin(),w.end(),compare);
00093                 for(int j=0;j<n;j++){
00094                         int k = w[ny - j -1].second;
00095                         weights[i][k] = true ;
00096                 }
00097                 
00098         }
00099         for(int i=0;i<ny;i++){
00100                 vector< pair<real,int> > w(nx);
00101                 for(int j=0;j<nx;j++){
00102                         w[j].first=nkmat[j][i];
00103                         w[j].second=j;
00104                 }
00105                 sort(w.begin(),w.end(),compare);
00106                 for(int j=0;j<n;j++){
00107                         int k = w[nx - j -1].second;
00108                         weights[k][i] = true ;
00109                 }
00110                 
00111         }
00112         
00113         
00114 }
00115 
00119 void static extractWeightList(const vector< vector<bool> >& wfilter , vector< pair<int,int> >& wlist){
00120         int nx = wfilter.size();
00121         int count = 0;
00122         for(int i=0;i<nx;i++){
00123                 int ny = wfilter[i].size();
00124                 
00125                 for(int j=0;j<ny;j++){
00126                         if(wfilter[i][j]==true){
00127                                 pair<int,int> w;
00128                                 w.first = i;
00129                                 w.second = j;
00130                                 wlist.push_back(w);
00131                                 count++;
00132                         }
00133                 }
00134         }
00135         //cout<<count<<" entries added to wlist"<<endl;
00136 }
00137 
00141 void static calculateEuclDist(const Mat& coords,Mat& dist){
00142         int n = coords.nrows();
00143         dist = Mat(n,n,0.0);
00144         for(int i=0;i<n;i++){
00145                 for(int j=0;j<n;j++){
00146                         for(int k=0;k<3;k++){
00147                                 dist[i][j] += (coords[i][k]-coords[j][k])*(coords[i][k]-coords[j][k]);
00148                         }
00149                         dist[i][j] = sqrt(dist[i][j]);
00150                 }
00151         }       
00152 }
00153 
00154 
00155 
00156 
00167 real static calcTransformation4(const Mat &xmat,const Mat& ymat,const Mat& wij,const Mat &nk,Mat& rot,Vec& xm,Vec& ym){
00168         int newn = xmat.nrows()+ymat.nrows();
00169         Mat xmat2(newn,xmat.ncols());
00170         Mat ymat2(newn,xmat.ncols());
00171         Vec weights(newn);
00172         for(int i=0;i<xmat.nrows();i++){
00173                 int max=0;
00174                 for(int j=0;j<ymat.nrows();j++){
00175                         if(wij[i][max]<wij[i][j]) max=j;
00176                 }
00177                 weights[i]=wij[i][max];
00178                 for(int j=0;j<3;j++){
00179                         xmat2[i][j]=xmat[i][j];
00180                         ymat2[i][j]=ymat[max][j];
00181                 }
00182         }
00183         for(int i=0;i<ymat.nrows();i++){
00184                 int max=0;
00185                 for(int j=0;j<xmat.nrows();j++){
00186                         if(wij[max][i]<wij[j][i]) max=j;
00187                 }
00188                 weights[i+xmat.nrows()]=wij[max][i];
00189                 for(int j=0;j<3;j++){
00190                         xmat2[i+xmat.nrows()][j]=xmat[max][j];
00191                         ymat2[i+xmat.nrows()][j]=ymat[i][j];
00192                 }
00193         }
00194         real sum = 0.0;
00195         for(int i=0;i<xmat.nrows()+ymat.nrows();i++){
00196                 sum += weights[i];
00197         }
00198         real error;
00199         rot = Mat(3,3);
00200         xm = weightedCentroid(xmat,weights);
00201         ym = weightedCentroid(ymat,weights);
00202         xmat2 -= xm;
00203         ymat2 -=ym;
00204         rot = weightedRotationFromMatchedPoints(xmat2,ymat2,weights,error);
00205         error /= sum;
00206         cout<<"computation complete "<<error<<endl;
00207         return error;
00208 }
00209 
00210 
00211 
00212 
00213 
00214 
00218 void static autoThreshLP(const Mat& dist1,const Mat& dist2,const Mat& nk,const vector< pair<int,int> >& wlist,const vector< vector<bool> >& wfilter,Mat& wm){
00219         int nterms = 0;
00220         int n = wlist.size();
00221         int rows = dist1.nrows();
00222         int cols = dist2.nrows();
00223         int *ia = new int[1+50000];
00224         int *ja = new int[1+50000];
00225         double *ar = new double[1+50000];
00226         const real sigma = 0.07;
00227         vector< pair<int,int> > bpairs;
00228         vector<real> products;
00229         const double threshs[] = {0.7,0.75,0.8,0.85,0.9,0.95,0.97,0.98,0.99,0.995,0.999};
00230         const int threshs_size = 11;
00231         real thresh = threshs[0];
00232         for(int i=0;i<n;i++){
00233                 for(int j=i+1;j<n;j++){
00234                         int xa = wlist[i].first;
00235                         int ya = wlist[i].second;
00236                         int xb = wlist[j].first;
00237                         int yb = wlist[j].second;
00238                         if(nk[xa][ya]>thresh && nk[xb][yb]>thresh){
00239                                 real d1 = dist1[xa][xb];
00240                                 real d2 = dist2[ya][yb];
00241                                 real diff = (d1-d2)*(d1-d2);
00242                                 diff = diff/sigma*sigma;
00243                                 real ekernel = exp(-diff);
00244                                 real product = ekernel*nk[xa][ya]*nk[xb][yb];
00245                                 if(product>thresh){
00246                                         nterms++;
00247                                         bpairs.push_back(pair<int,int>(i,j));
00248                                         products.push_back(product);
00249                                 }
00250                         }       
00251                 }
00252         }
00253         //cout<<"number of quadratic constraints "<<nterms<<endl;
00254         for(int i=1;i<threshs_size;i++){
00255                 if(nterms<0.15*n) break;
00256                 thresh = threshs[i];
00257                 nterms = 0;
00258                 for(unsigned int j=0;j<products.size();j++){
00259                         if(products[j]>thresh) nterms++;
00260                 }       
00261         }
00262 //      cout<<"threshold selected "<<thresh<<endl;
00263 //      cout<<"terms to be inserted "<<nterms<<endl;
00264         int count = 1;
00265         vector<real> row_bnds;
00266         for(unsigned int k=0;k<products.size();k++){
00267                 if(products[k]>thresh){
00268                         int i = bpairs[k].first;
00269                         int j = bpairs[k].second;
00270                         int rownumber = rows + cols + (count+1)/2;
00271             if (2 * n + count >= 50000) cout << "overflow"  << endl; 
00272                         ia[2*n + count] = rownumber;
00273                         ja[2*n + count] = i+1;
00274                         ar[2*n + count] = 1.0;
00275                         count++;
00276                         ia[2*n + count] = rownumber;
00277                         ja[2*n + count] = j+1;
00278                         ar[2*n + count] = 1.0;  
00279                         row_bnds.push_back(1+products[k]);
00280                         count++;
00281                 }
00282         }
00283         //nconst = number of constraints
00284         int nconst = rows+cols + nterms;
00285         LPX *lp;
00286         lp = lpx_create_prob();
00287         lpx_set_int_parm(lp,LPX_K_MSGLEV,1);
00288         lpx_set_obj_dir(lp,LPX_MIN);
00289         lpx_add_cols(lp,n);
00290         lpx_add_rows(lp,nconst);
00291         for(int i=0;i<n;i++){
00292                 lpx_set_col_bnds(lp,i+1,LPX_DB,0.0,1.0);
00293         }
00294         //each weight appears in 2 constraints
00295         //sconst : size of constraint vector
00296         int sconst = 2*n + 2*nterms;    
00297                 
00298         //generate row contraints
00299         count = 1;
00300         for(int i=0;i<rows;i++){
00301                 for(int j=0;j<cols;j++){
00302                         if(wfilter[i][j]){
00303                                 ia[count] = i+1;
00304                                 ja[count] = count;
00305                                 ar[count] = 1.0;
00306                                 ia[n+count] = j+rows+1;
00307                                 ja[n+count] = count;
00308                                 ar[n+count] = 1.0;
00309                                 count++;        
00310                         }       
00311                 }
00312         }
00313         //generate the quadratic constraints
00314         for(int i=0;i<nterms;i++){
00315                 lpx_set_row_bnds(lp,rows+cols+i+1,LPX_LO,row_bnds[i],0.0);
00316         }
00317         
00318         for(int i=0;i<rows+cols;i++){
00319                 lpx_set_row_bnds(lp,i+1,LPX_LO,1.0,0.0);
00320         }
00321         lpx_load_matrix(lp,sconst,ia,ja,ar);
00322         //generate the linear coefficients of the weights which is just the nodekernels
00323         for(int i=0;i<n;i++){
00324                 int a = wlist[i].first;
00325                 int b = wlist[i].second;
00326                 real coef = 1.0 - nk[a][b]; //+ lambda*abs(xcdist[a]-ycdist[b]);
00327                 lpx_set_obj_coef(lp,i+1,coef);
00328         }
00329 //      time_t t1 = time(NULL);
00330         lpx_simplex(lp);
00331 //      time_t t2 = time(NULL);
00332 //      cout<<difftime(t2,t1)<<endl;
00333         
00334         wm = Mat(rows,cols,0.0);
00335         for(int i=0;i<n;i++){
00336                 wm[wlist[i].first][wlist[i].second]=lpx_get_col_prim(lp,i+1);
00337         }
00338         lpx_delete_prob(lp);
00339         free(ia);
00340         free(ja);
00341         free(ar);               
00342 }
00343 
00351 void static performLP(PMolecule name1,MoleculeTemplate name2,Mat& wm,bool isweighted){
00352 try{    
00353 //      cout<<"performing lp on "<<name1<<" "<<name2<<endl;
00354         Mat xprpt,yprpt;
00355 
00356     xprpt = name1->chem ; 
00357     yprpt = name2->chem ; 
00358     
00359         Mat dev;
00360 //      if(isweighted){
00361 //              load(name2+"Dev.mat",dev);
00362 //      }else{
00363                 dev=Mat(yprpt.nrows(),yprpt.ncols(),1.0);
00364 //      }
00365         vector<int> common_prpt;
00366         for(int i=0;i<yprpt.ncols();i++){
00367                 if(yprpt[0][i]<5000.0 && xprpt[0][i]<5000.0){
00368                         common_prpt.push_back(i);
00369                 }       
00370         }
00371         Mat xprpt2(xprpt.nrows(),common_prpt.size(),0.0);
00372         Mat yprpt2(yprpt.nrows(),common_prpt.size(),0.0);
00373         Mat dev2(yprpt.nrows(),common_prpt.size(),0.0);
00374         for(unsigned int j=0;j<common_prpt.size();j++){
00375                 int k = common_prpt[j];
00376                 for(int i=0;i<xprpt.nrows();i++){
00377                         xprpt2[i][j]=xprpt[i][k];
00378                 }
00379                 for(int i=0;i<yprpt.nrows();i++){
00380                         yprpt2[i][j]=yprpt[i][k];
00381                         dev2[i][j]=dev[i][k];
00382                 }
00383         }
00384         Mat nk;
00385         nodekernel(xprpt2,yprpt2,dev2,0.3*common_prpt.size(),nk);
00386         const int rows=nk.nrows();
00387         const int cols=nk.ncols();
00388         const int yjperxi = 5;
00389         vector< vector<bool> > wfilter;
00390         vector< pair<int,int> > wlist;
00391         findRelevantWeights(nk,wfilter,yjperxi);
00392         extractWeightList(wfilter,wlist);
00393         const int n = wlist.size();     
00394         
00395         /*reading distances*/
00396         Mat dist1,dist2;
00397         calculateEuclDist(name1->geom,dist1);
00398         calculateEuclDist(name2->geom,dist2);
00399         if(dist1.nrows()!=rows || dist2.nrows()!=cols){
00400                 cout<<dist1.nrows()<<" "<<dist2.nrows()<<" "<<rows<<" "<<cols<<endl;
00401                 PLERROR("the dimensions in dist files and nkmat do not match\n");
00402         }
00403         Mat xmat,ymat;
00404 
00405     xmat = name1->geom ; 
00406     ymat = name2->geom ; 
00407     
00408         //cout<<"got vertex coords"<<endl;
00409         for(int i=0;i<n;i++){
00410                 int a = wlist[i].first;
00411                 int b = wlist[i].second;
00412                 if(a>=rows || b>=cols){
00413                         cout<<i<<" "<<a<<" "<<b<<endl;
00414                 }
00415         }
00416         //calcLinearWeights(dist1,dist2,nk,wlist,wfilter,sigma,thresh,wm);
00417         autoThreshLP(dist1,dist2,nk,wlist,wfilter,wm);
00418         //calcTransformation4(xmat,ymat,wm,nk,rot,xm,ym);
00419         //writeAlignment(name1,xm,ym,rot,"lpw4"); 
00420 }catch(PLearnError e){
00421         cout<<e.message()<<endl;
00422 }       
00423 }
00424 
00425 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines