PLearn 0.1
|
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