From 7bc94eb6c7d86dfaddea5ba563cede99d0c5cf98 Mon Sep 17 00:00:00 2001 From: giopina Date: Mon, 4 Jul 2016 16:14:15 +0200 Subject: [PATCH] first version --- Bead.h | 220 +++++ BeadList.h | 115 +++ Connection.h | 59 ++ CovMatrix1d.cc | 347 ++++++++ ElasticNet.cc | 1168 ++++++++++++++++++++++++++ ElasticNet.h | 160 ++++ EntrComp.C | 138 +++ Makefile | 97 +++ Matrix.cc | 1125 +++++++++++++++++++++++++ Matrix.h | 217 +++++ NormalModes.cpp | 111 +++ NormalModes.h | 39 + PrincipalComp.cc | 671 +++++++++++++++ PrincipalComp.h | 99 +++ RiboGM.C | 60 ++ Structure3d.cpp | 649 ++++++++++++++ Structure3d.h | 110 +++ Vector3d.h | 129 +++ ang.C | 499 +++++++++++ ang_C1p_P.C | 398 +++++++++ common.h | 44 + compdist.C | 167 ++++ correlation.cc | 46 + countn.C | 87 ++ distC2.C | 534 ++++++++++++ enm_new.C | 464 ++++++++++ fast_Bhatt.C | 187 +++++ fast_rwsip.C | 147 ++++ fit.C | 242 ++++++ genConf.C | 175 ++++ io.h | 81 ++ kabsch.cc | 530 ++++++++++++ kabsch.h | 19 + lapack_matrix_routines_wrapper_v3.cc | 252 ++++++ my_malloc.h | 870 +++++++++++++++++++ myjacobi.cc | 105 +++ pca.C | 251 ++++++ qcprot.cpp | 455 ++++++++++ qcprot.h | 164 ++++ quaternion.c++ | 342 ++++++++ quaternion.h | 253 ++++++ trace.C | 37 + traj2com.C | 54 ++ vectop.cc | 187 +++++ 44 files changed, 12104 insertions(+) create mode 100644 Bead.h create mode 100644 BeadList.h create mode 100644 Connection.h create mode 100644 CovMatrix1d.cc create mode 100644 ElasticNet.cc create mode 100644 ElasticNet.h create mode 100644 EntrComp.C create mode 100644 Makefile create mode 100644 Matrix.cc create mode 100644 Matrix.h create mode 100644 NormalModes.cpp create mode 100644 NormalModes.h create mode 100644 PrincipalComp.cc create mode 100644 PrincipalComp.h create mode 100644 RiboGM.C create mode 100644 Structure3d.cpp create mode 100644 Structure3d.h create mode 100644 Vector3d.h create mode 100644 ang.C create mode 100644 ang_C1p_P.C create mode 100644 common.h create mode 100644 compdist.C create mode 100644 correlation.cc create mode 100644 countn.C create mode 100644 distC2.C create mode 100644 enm_new.C create mode 100644 fast_Bhatt.C create mode 100644 fast_rwsip.C create mode 100644 fit.C create mode 100644 genConf.C create mode 100644 io.h create mode 100644 kabsch.cc create mode 100644 kabsch.h create mode 100644 lapack_matrix_routines_wrapper_v3.cc create mode 100644 my_malloc.h create mode 100644 myjacobi.cc create mode 100644 pca.C create mode 100644 qcprot.cpp create mode 100644 qcprot.h create mode 100644 quaternion.c++ create mode 100644 quaternion.h create mode 100644 trace.C create mode 100644 traj2com.C create mode 100644 vectop.cc diff --git a/Bead.h b/Bead.h new file mode 100644 index 0000000..003045a --- /dev/null +++ b/Bead.h @@ -0,0 +1,220 @@ +/* + * Bead.h + * + * Created on: Jul 19, 2012 + * Author: gpolles + * + * Bead class. + * Each bead has 3 coordinates and pointers to its neighbors + * and to a domain. + * It have also some information from the pdb file + * as chainId, and betaFactor + * + */ + +#ifndef BEAD_H_ +#define BEAD_H_ + +#include + +#include "Vector3d.h" +#include "common.h" + + +class RigidDomain; + +class Bead +{ +public: + Bead() : _index(0), _num_neighbors(0),_num_covBonded(0), _domain(NULL), + _isEdge(false), _coordinates(Vector3d(0,0,0)), + _atomType(""),_resNum(0),_resName(""),_hetatm(false), + _betaFactor(0.0), _pdbPosition(0), _chainId(0) {} + virtual ~Bead(){} + + double distance(Bead* b){ + return _coordinates.distance(b->_coordinates); + } + + double distance(Bead& b){ + return _coordinates.distance(b._coordinates); + } + + double distanceSQ(Bead* b){ + return _coordinates.distanceSQ(b->_coordinates); + } + + double distanceSQ(Bead& b){ + return _coordinates.distanceSQ(b._coordinates); + } + + RigidDomain* getDomain() const { + return _domain; + } + + void setDomain(RigidDomain* domain){ + _domain = domain; + } + + int getIndex() const { + return _index; + } + + bool getHet() { + return _hetatm; + } + + void setHet(bool HET) { + _hetatm=HET; + } + + void setIndex(int index) { + _index = index; + } + + std::vector& getNeighbors() { + return _neighbors; + } + + void setNeighbors(std::vector neighbors) { + _neighbors = neighbors; + } + + int getNumNeighbors() const { + return _num_neighbors; + } + + void setNumNeighbors(int numNeighbors) { + _num_neighbors = numNeighbors; + } + + void addNeighbor(Bead* neighbor){ + _neighbors.push_back(neighbor); + _num_neighbors = _neighbors.size(); + } + + + + /************************/ +std::vector& getCovBonded() { + return _covBonded; + } + + void setCovBonded(std::vector covBonded) { + _covBonded = covBonded; + } + + int getNumCovBonded() const { + return _num_covBonded; + } + + void setNumCovBonded(int numCovBonded) { + _num_covBonded = numCovBonded; + } + + void addCovBond(Bead* covbond){ + _covBonded.push_back(covbond); + _num_covBonded = _covBonded.size(); + } + + /***********************/ + + + bool isEdge(){ + return _isEdge; + } + + void setEdge(bool b){ + _isEdge =b; + } + + double getBetaFactor() const { + return _betaFactor; + } + + void setBetaFactor(double betaFactor) { + _betaFactor = betaFactor; + } + + //##################################### + std::string getAtomType() { + return _atomType; + } + void setAtomType(std::string type) { + _atomType=type; + } + + int getBeadType() { + return _beadType; + } + void setBeadType(int type) { + _beadType=type; + } + + std::string getResName() { + return _resName; + } + void setResName(std::string name) { + _resName=name; + } + + int getResNum() { + return _resNum; + } + void setResNum(int resnum) { + _resNum=resnum; + } + //#################################### + + + + + const Vector3d& getCoordinates() const { + return _coordinates; + } + + void setCoordinates(const Vector3d& coordinates) { + _coordinates = coordinates; + } + + int getPdbPosition() const { + return _pdbPosition; + } + + void setPdbPosition(int pdbPosition) { + _pdbPosition = pdbPosition; + } + + int getChainId() const { + return _chainId; + } + + void setChainId(int chainId) { + _chainId = chainId; + } + +private: + int _index; + int _num_neighbors; + int _num_covBonded; + RigidDomain* _domain; + std::vector _neighbors; + std::vector _covBonded; + bool _isEdge; + Vector3d _coordinates; + double _betaFactor; + int _pdbPosition; + bool _hetatm; + + + int _beadType; // 0 -> P + // 1 -> sugar + // 2 -> base + std::string _atomType; + std::string _resName; + int _resNum; + + int _chainId; +}; + +#endif /* BEAD_H_ */ diff --git a/BeadList.h b/BeadList.h new file mode 100644 index 0000000..1974a67 --- /dev/null +++ b/BeadList.h @@ -0,0 +1,115 @@ +/* + * BeadList.h + * + * Created on Mar 25, 2014 + * Author: Giovanni Pinamonti + * + * BeadList class. + * + */ + +#ifndef BEADLIST_H_ +#define BEADLIST_H_ + +#include +#include "Bead.h" +#include + +class BeadList +{ + + private: + std::vector _atomTypes; + std::vector _resNumbers; + bool _noH; + bool _n_1_9; + + public: + //creatori + BeadList():_noH(true), _n_1_9(false) {}; + + BeadList(std::vector namelist):_noH(true), _n_1_9(false) { + if(namelist.size()>0){ + for (int i=0;i numlist):_noH(true), _n_1_9(false) { + _resNumbers=numlist; + } + /* + //funzioni + */ + void setResNum(std::vector numlist){ + // cout<<"cacca"< types){ + // std::vector + _atomTypes=types; + } + + + inline bool check(Bead b){ + + // cout<<"cacca"<0)||(_n_1_9)){ + bool lacacca=false; + std::string s1=b.getAtomType(); + for (int i=0;i<_atomTypes.size();++i){ + std::string s2=_atomTypes.at(i); + std::replace( s2.begin(), s2.end(), 'p', '\''); + // cout<0){ + bool lacacca=false; + int r=b.getResNum(); + for (int i=0;i<_resNumbers.size();++i){ + if (r==_resNumbers.at(i)) {lacacca=true;} + } + if(!lacacca) return false; + } + + // cout<<_noH< + +class Connection +{ + public: + Connection(int i, int j) { + _i1=i; + _i2=j; + //QUI DOVREI CONTROLLARE CHE SIANO POSITIVI E DIVERSI + } + + inline int GetI1() {return _i1;} + inline int GetI2() {return _i2;} + + inline bool SetI1(int i) { if((i>0)&&(i!=_i2)){ _i1=i; return true;} return false; /*CONTROLLA CHE SIANO POSITIVI E DIVERSI*/} + inline bool SetI2(int i) { if((i>0)&&(i!=_i1)){ _i2=i; return true;} return false; /*CONTROLLA CHE SIANO POSITIVI E DIVERSI*/} + + private: + int _i1; + int _i2; +}; + + +//######################################################## +//### funzioni utili che utilizzano la suddetta classe ### +//######################################################## + +inline bool RandSwap(Connection &p,Connection &q) { + double r=((double) rand() / ((double) RAND_MAX)); + if(r>0.5){ + int i1=p.GetI1(); + int j1=q.GetI1(); + if( (p.SetI1(j1)) && + (q.SetI1(i1)) ) + return true; + } + else{ + // int j2=q.GetI2(); + int j1=q.GetI1(); + int i2=p.GetI2(); + if( (p.SetI2(j1)) && + (q.SetI1(i2)) ) + return true; + } + return false; +} + +#endif /* CONNECTION_H_ */ diff --git a/CovMatrix1d.cc b/CovMatrix1d.cc new file mode 100644 index 0000000..1b49029 --- /dev/null +++ b/CovMatrix1d.cc @@ -0,0 +1,347 @@ +/* + Created by Giovanni Pinamonti + PhD student @ + Scuola Internazionale Superiore di Studi Avanzati, Trieste, Italy + November 21st 2013 + */ + +//extern "C" { +//#include +//} + +#define TOL 0.000001 + +#include +#include +#include + +#include "Matrix.h" +#include "my_malloc.h" +#include "io.h" +//#include "Vector3d.h" + +#include "lapack_matrix_routines_wrapper_v3.cc" + +using namespace std; + +//Constructors for the class CovMatrix1d +CovMatrix1d:: CovMatrix1d():CovMatrix(){ //forse sta cosa non si dovrebbe fare. Se non dichiaro la dimensione della matrice non la creo... + _DIM=1; +}//END OF CONSTRUCTOR + +CovMatrix1d:: CovMatrix1d(int N):CovMatrix(){ + //size of matrix + _DIM=1; + _N=N; + _size=_DIM*_N; + _matrix=d2t(_size,_size); +}//END OF CONSTRUCTOR + + +CovMatrix1d:: CovMatrix1d(int N, double **matrix):CovMatrix(){ + _DIM=1; + _N=N; + _size=_DIM*_N; + _matrix=d2t(_size,_size); + for(int i=0; i<_size; ++i){ + for(int j=0; j<_size; ++j){ + _matrix[i][j]=matrix[i][j]; + }//enddo + }//enddo +}//END OF CONSTRUCTOR + + +//############################################# + +void CovMatrix1d::SetEigenvectors(std::vector > evec){ + int numModes=evec.size(); + int dim=evec.at(0).size(); + if(_size!=numModes){ + cerr<0) _eigenvectors.clear(); + + try{ + _eigenvectors.resize(numModes); + for (int i = 0; i < numModes; ++i) { + _eigenvectors[i].reserve(dim); + for (int j = 0; j < dim; ++j) { + double p1d=evec.at(i).at(j); + _eigenvectors[i].push_back(p1d); + } + } + + }catch(std::bad_alloc const&){ + cerr << endl << "NormalMode eigenvector memory allocation failed!" << endl << flush; + exit(1); + } +} + + +void CovMatrix1d::Compose(){ + if(_size>0) free_d2t(_matrix);// AM I SURE THAT THIS IS SAFE?? + _size=_eigenvalues.size(); + if(_size<1){ + cout<<"Matrix.Compose --> ERROR: unvalid size of the matrix. Maybe you forgot to set the eigenvalues"< ERROR: unvalid size of the matrix. Maybe you forgot to set the eigenvectors"<= k 0 no */ + for(int k = 0;k< _size; k++){ + // imnotdoingthis?? /* Run the sum backwards to minimize roundoff errors in the sum of inverse eigenvalues */ + // if (fabs(_eigenvalues[k]) < tol) continue; + + _matrix[i][j] += _eigenvalues[k]*_eigenvectors.at(k).at(i)*_eigenvectors.at(k).at(j); + + } + if (i!=j){ /* symmetry */ + _matrix[j][i]=_matrix[i][j]; + }//endif + }//enddo j + }//endoi +}//end function + +void CovMatrix1d::Compose(std::vector eval, std::vector > evec){ + SetEigenvalues(eval); + SetEigenvectors(evec); + _swapEigen(); + return Compose(); +} + +//############################################# + +void CovMatrix1d:: Decompose(){ + + if(_eigenvalues.size()>0){ + _eigenvalues.clear(); + _eigenvectors.clear(); + cout<<"deleting old eigenstuff"< eigenvec; + for(int j=0;j<_N;++j){ + double ivec=temp_eigenvectors[i][j]; + eigenvec.push_back(ivec); + } + _eigenvalues.push_back(temp_eigenvalues[i]); + _eigenvectors.push_back(eigenvec); + } + free_d1t(temp_eigenvalues); + free_d2t(temp_eigenvectors); +} + +//TODO you have to decompose before invert. Can I modify this so that it check if it was done? +Matrix CovMatrix1d::GetInverse(){ + //TODO: sta cosa non so se mi serve in fondo... + Matrix tempMat(0); + return tempMat; + // double **invmat=d2t(_size,_size); + // double *temp_eigenvalues=d1t(_size); + // double **temp_eigenvectors=d2t(_size,_size); + // for(int i=0;i<_size;++i){ + // temp_eigenvalues[i]=_eigenvalues.at(i); + // for(int j=0;j<_N;++j){ + // temp_eigenvectors[i][_DIM*j+0]=_eigenvectors.at(i).at(j).X; + // temp_eigenvectors[i][_DIM*j+1]=_eigenvectors.at(i).at(j).Y; + // temp_eigenvectors[i][_DIM*j+2]=_eigenvectors.at(i).at(j).Z; + // } + // } + // spectral_singular_inversion(invmat,_size,temp_eigenvalues,temp_eigenvectors,TOL); + // free_d1t(temp_eigenvalues); + // free_d2t(temp_eigenvectors); + // // cout<<"caccacca"<= _size) break; + sprintf(filename,"%s_eigenvector_%d.dat",name,i); + fp =open_file_w(filename); + for(int j=0; j < _N; j++){ + fprintf(fp,"%4d ",i); + fprintf(fp,"%lf ",_eigenvectors.at(_size-i-1).at(j)); + fprintf(fp,"\n"); + } + fclose(fp); + printf("Eigenvector number %4d written to file %s\n",i,filename); + } +} + + +void CovMatrix1d::dumpMatrix(const char *name){ + char filename[200]; + // FILE *fp; + sprintf(filename,"%s_matrix.dat",name); + ofstream fp; + fp.open(filename); + for(int i=0; i < _N; i++){ + for(int j=0; j < _N; j++){ + fp< tempvec; + + for(int i=0; i<_size/2; ++i){ + // switch the eigenvalues + tempval=_eigenvalues[i]; + _eigenvalues[i]=_eigenvalues[_size-1-i]; + _eigenvalues[_size-1-i]=tempval; + //switch the eigenvectors + tempvec=_eigenvectors[i]; + _eigenvectors[i]=_eigenvectors[_size-1-i]; + _eigenvectors[_size-1-i]=tempvec; + } +} + + + +//#################################################################### +//#################################################################### +//#################################################################### + +void CovMatrix1d::dumpMSF(const char*name){ + char filename[200]; + FILE *fp; + sprintf(filename,"%s_mean_square_displ.dat",name); + fp =open_file_w(filename); + for(int i=0; i < _N; i++){ + double temp=GetMSF(i); + fprintf(fp,"%4d %e\n",i,temp); + } + fclose(fp); + printf("Beads' mean square displacement written to file %s\n",filename); +} + + +double CovMatrix1d::GetMSF(int i){ + double temp= _matrix[i][i]; + return temp;} + + +//TODO: non so se mi serve sta cosa +//IntMatrix CovMatrix::GetInverse(){ + // double **invmat=d2t(_size,_size); + // double *temp_eigenvalues=d1t(_size); + // double **temp_eigenvectors=d2t(_size,_size); + // for(int i=0;i<_size;++i){ + // temp_eigenvalues[i]=_eigenvalues.at(i); + // for(int j=0;j<_N;++j){ + // temp_eigenvectors[i][_DIM*j+0]=_eigenvectors.at(i).at(j).X; + // temp_eigenvectors[i][_DIM*j+1]=_eigenvectors.at(i).at(j).Y; + // temp_eigenvectors[i][_DIM*j+2]=_eigenvectors.at(i).at(j).Z; + // } + // } + // spectral_singular_inversion(invmat,_size,temp_eigenvalues,temp_eigenvectors,TOL); + // free_d1t(temp_eigenvalues); + // free_d2t(temp_eigenvectors); + // IntMatrix tempMat(_N,invmat); + // free_d2t(invmat); + // return tempMat; +//} + + +void CovMatrix1d::Reduce(int NTOP, const CovMatrix1d&D){ + // reduces the matrix using the first NTOP eigenvectors of matrix D; + // Matrix D has to be already decomposed; + // if D.GetEigenval + double** temp_mat=d2t(NTOP,NTOP); + + for(int j=0;j CQ_m_j; + double * CQ_M_j=d1t(_size); + for(int m=0;m<_size;++m){ + for(int l=0;l<_N;++l){ + CQ_M_j[m]+=_matrix[m][l]*D._eigenvectors[_size-j-1].at(l); + }//enddo l + }//enddo m + for(int i=0;i0){ + _eigenvectors.clear(); + _eigenvalues.clear(); + } + for(int i=0;i<_size;++i){ + for(int j=0;j<_size;++j){ + _matrix[i][j]=temp_mat[i][j]; + } + } +} + + +double RWSIP(const CovMatrix1d &A,const CovMatrix1d &B){ + double RWSIP=0.0; + cout<<"### RWSIP ###"< evec_A=A._eigenvectors[m]; + for (int n = 6; n < nmodes; ++n) { + double eval_B=B._eigenvalues[n]; + vector evec_B=B._eigenvectors[n]; + double scal_prod=0; + for(int i=0; i +#include "omp.h" +#include "ElasticNet.h" +#include "io.h" +#include +#include +#include +#include "lapack_matrix_routines_wrapper_v3.cc" +#include "BeadList.h" + +#define BROKEN_CHAIN_CUTOFF 7.0 + +#define DEFAULT_CUTOFF 10.0 +#define DEFAULT_NTOP 10 + +//###### COSTRUTTORI DELLA CLASSE ###### +ElasticNet::ElasticNet(): + // _covMatrix(0) + _N(0), + _ntop_vectors(DEFAULT_NTOP), + _cutoff(DEFAULT_CUTOFF), + + _R_BB(-1), + _R_B(-1), + _R_PP(-1), + _R_prot(-1), + _R_N7(-1), + + _TRACCIAMENTO(false), + _Kcov(1.0), + _KBB(1.0), + _KB(1.0), + _KPP(1.0), + _COM(false), + _OUTRES(false), + // _PROTEIN(false), + _FAST(false), + _DPF(false), + _RANDOM(false), + _RARAND(false), + _EXP(false), + _DUMPMODES(false), + _DUMPEIGENVALUES(true), + _DUMPCOVMAT(false), + _DUMPDISTFLUC(false), + _DUMPREDCOVMAT(false), + _DUMPMSD(true) +{ + _COM=false; +} + +ElasticNet::ElasticNet(const char *filename): + _ntop_vectors(DEFAULT_NTOP), + _cutoff(DEFAULT_CUTOFF), + + _R_BB(-1), + _R_B(-1), + _R_PP(-1), + _R_prot(-1), + _R_N7(-1), + + _TRACCIAMENTO(false), + _Kcov(1.0), + _KBB(1.0), + _KB(1.0), + _KPP(1.0), + _COM(false), + _OUTRES(false), + // _PROTEIN(false), + _FAST(false), + _DPF(false), + _RANDOM(false), + _RARAND(false), + _EXP(false), + _DUMPMODES(false), + _DUMPEIGENVALUES(true), + _DUMPCOVMAT(false), + _DUMPDISTFLUC(false), + _DUMPREDCOVMAT(false), + _DUMPMSD(true) +{ + // cout<<"ElasticNet::reading file "<Decompose() + if(_TRACCIAMENTO) { + if(tracciamento()<0){cout<<"ERROR IN TRACCIAMENTO"< decompose..."< _covMatrix.dumpEigenvalues()"< _covMatrix.dumpEigenvectors()"< Structure"< dumpMSF()()"< PROBLEM: won't dump principal modes in FAST configuration!"< contatto normale + // 3 -> legame covalente + // 10 -> contatto proteina-proteina + // 11 -> contatto proteina-RNA + + // double R_BB=5.0; + + if(_R_B<0) _R_B =_cutoff; + if(_R_PP<0) _R_PP =_cutoff; + if(_R_BB<0) _R_BB =_cutoff; + if(_R_prot<0) _R_prot=_cutoff; + if(_R_N7<0) _R_N7 =_cutoff; + + cout<<" *** RADIUS OF CUT-OFF *** "< base-base + //20 -> base other + //99 -> P-P + //77 -> prot-prot + //77 -> prot-RNA + + + if (bi.getAtomType()=="N7"){ + // if (bj.getAtomType()=="N7"){ + // Rc=_R_BB; kind_o_int=22;} + // else + {Rc=_R_N7; kind_o_int=20;} + } + if (bj.getAtomType()=="N7"){ + // if (bi.getAtomType()=="N7"){ + // Rc=_R_BB; kind_o_int=22;} + // else + {Rc=_R_N7; kind_o_int=20;} + } + + + // if ((bi.getAtomType()=="C2")||(bi.getAtomType()=="N7")){ + // if ((bj.getAtomType()=="C2")||(bj.getAtomType()=="N7")){ + if (bi.getAtomType()=="C2"){ + if (bj.getAtomType()=="C2"){ + Rc=_R_BB; kind_o_int=22;} + else + {Rc=_R_B; kind_o_int=20;} + } + // if ((bj.getAtomType()=="C2")||(bj.getAtomType()=="N7")){ + // if ((bi.getAtomType()=="C2")||(bi.getAtomType()=="N7")){ + if (bj.getAtomType()=="C2"){ + if (bi.getAtomType()=="C2"){ + Rc=_R_BB; kind_o_int=22;} + else + {Rc=_R_B; kind_o_int=20;} + } + + + + if (bi.getAtomType()=="P") + if (bj.getAtomType()=="P"){ + Rc=_R_PP; kind_o_int=99;} + + double Rhybrid=(_cutoff+_R_prot)*0.5; + //protein RNA interaction + if (bi.getAtomType()=="CA"){ + if (bj.getAtomType()=="CA"){ + Rc=_R_prot;kind_o_int=77;} + else{ Rc=Rhybrid;kind_o_int=70;} + } + else{ + if (bj.getAtomType()=="CA"){ + Rc=Rhybrid;kind_o_int=70;} + } + + + if (d < Rc) _cmap[i][j]=_cmap[j][i]=kind_o_int; + }//end do j + }//end do i + + + for(int i=0; i < _N-1; i++){ + /* ### stiffest springs between bonded nodes ### */ + Bead bi=_structure.getBead(i); + if (bi.getAtomType()=="P"){ + int j=i+1; + Bead bj=_structure.getBead(j); + if (bj.getAtomType()=="C1'"){ + double d=_structure.getDistance(i,j); + if (d < BROKEN_CHAIN_CUTOFF) + _cmap[i][j]=_cmap[j][i]=3; + }//endif + }//endif P + + else if (bi.getAtomType()=="C1'"){ + int j=i+1; + Bead bj=_structure.getBead(j); + if (bj.getAtomType()=="C2"){ + double d=_structure.getDistance(i,j); + if (d < BROKEN_CHAIN_CUTOFF) + _cmap[i][j]=_cmap[j][i]=3; + }//endif + if (bj.getAtomType()=="N7"){ + double d=_structure.getDistance(i,j); + if (d < BROKEN_CHAIN_CUTOFF) + _cmap[i][j]=_cmap[j][i]=3; + }//endif + + j=i+2; + if(j>_N-1) continue; + bj=_structure.getBead(j); + if (bj.getAtomType()=="P"){ + double d=_structure.getDistance(i,j); + if (d < BROKEN_CHAIN_CUTOFF) + _cmap[i][j]=_cmap[j][i]=3; + }//endif + + j=i+3; + if(j>_N-1) continue; + bj=_structure.getBead(j); + if (bj.getAtomType()=="P"){ + double d=_structure.getDistance(i,j); + if (d < BROKEN_CHAIN_CUTOFF) + _cmap[i][j]=_cmap[j][i]=3; + }//endif + + }//endif C1' + + else if (bi.getAtomType()=="N7"){ + int j=i+1; + Bead bj=_structure.getBead(j); + if (bj.getAtomType()=="C2"){ + double d=_structure.getDistance(i,j); + if (d < BROKEN_CHAIN_CUTOFF) + _cmap[i][j]=_cmap[j][i]=3; + }//endif + }//endif N7 + + }//enddo I + +}//end function + + +void ElasticNet::constructRandMat(double **intmat){ + cout<<"ElasticNet: creating a random interaction matrix"< connessioni; + for(int i=0;i<_N;++i){ + for (int j=i+1;j<_N;++j){ + if (i==j) continue; + if(_cmap[i][j]!=0){ + //n_conn++; + if(_cmap[i][j]>1) cout<0){cout<<"[i] were already connected"<0){cout<<"[j] were already connected"<n2swap) break; + }//end for + cout<<"Ho scambiato "<1) cout< Problem opening file: "<>stringa)) continue; + + if(!strncmp(stringa,"Int_Range",9)){iss>>argument; sscanf(argument,"%lf",&_cutoff);} + if(!strncmp(stringa,"R_BB",4)){iss>>argument; sscanf(argument,"%lf",&_R_BB);} + if(!strncmp(stringa,"R_Ba",4)){iss>>argument; sscanf(argument,"%lf",&_R_B);} + if(!strncmp(stringa,"R_PP",4)){iss>>argument; sscanf(argument,"%lf",&_R_PP);} + if(!strncmp(stringa,"R_prot",6)){iss>>argument; sscanf(argument,"%lf",&_R_prot);} + if(!strncmp(stringa,"R_N7",4)){iss>>argument; sscanf(argument,"%lf",&_R_N7);} + + + if(!strncmp(stringa,"N_SLOWEST_EIGENVECT",19)){iss>>argument; sscanf(argument,"%d",&_ntop_vectors);} + + + // if(!strncmp(stringa,"Smooth_Cutoff",13)) sscanf(argument,"%d",&Smooth_Cutoff); + // if(!strncmp(stringa,"nodesXnucl",10)) sscanf(argument,"%d",&nodesXnucl); + if(!strncmp(stringa,"K_COV",5)) {iss>>argument; sscanf(argument,"%lf",&(_Kcov)); cout<<"K_COV="<<_Kcov<>argument; sscanf(argument,"%lf",&(_KBB)); cout<<"K_BB="<<_KBB<>argument; sscanf(argument,"%lf",&(_KPP)); cout<<"K_PP="<<_KPP<>argument; sscanf(argument,"%lf",&(_KB)); cout<<"K_Ba="<<_KB<> name) + { + cout<> name) + { + cout<> num) + { + cout<>> Beads in the centers of mass for sugar and base <<< "<>> I'm gonna be faster than light! <<< "<>> Dump Partial Fluctuations <<< "<>> Dump covariance matrix <<< "<>> Dump distance fluctuations <<< "<>argument; sscanf(argument,"%lf",&seed); + srand(seed); + cout<<" @!#$$$&^%$! rAnDOm NeTWork &*#$&#%&^@## (seed = "<>argument; sscanf(argument,"%lf",&seed); + srand(seed); + cout<<" @!#$$$&^%$! rAnDOm NeTWork &*#$&#%&^@## (seed = "<>> exponential <<< "<>> Dump top modes <<< "<_beadList.size()) {cout<<"WARNING: OUTBEADS are more than beads!"< no need for that"< ERROR with the OUTBEADS list"<0){//Se ho calcolato covmatrix la uso ... + for(int mu=0;mu<3;mu++){ + for(int nu=0;nu<3;nu++){ + sigma+=d_ij[mu]*d_ij[nu]*( + _covMatrix.GetElement(3*i+mu,3*i+nu)+ + _covMatrix.GetElement(3*j+mu,3*j+nu)- + _covMatrix.GetElement(3*i+mu,3*j+nu)- + _covMatrix.GetElement(3*j+mu,3*i+nu)); + } + } + return sigma; + } + else{ // ... altrimenti ricalcolo solo gli elem che mi servono! + for(int mu=0;mu<3;mu++){ + for(int nu=0;nu<3;nu++){ + double C_ii_munu=0; + double C_jj_munu=0; + double C_ij_munu=0; + double C_ji_munu=0; + for(int alpha=0;alpha<3*_N-6;alpha++){ + double temp=_intMatrix.GetEigenval(alpha); + if (temp<0.000001) cout<<"WARNING: eigenvalue very small"< +#include "Connection.h" +#include "my_malloc.h" + +class ElasticNet +{ + + public: + ElasticNet(); + ElasticNet(Structure3d structure); + ElasticNet(const char *filename); //TODO costruttore direttamente con nome PDB + virtual ~ElasticNet(); + + CovMatrix & getCovMatrix(){ return (_covMatrix);} + IntMatrix & getIntMatrix(){ return (_intMatrix);} + + // double CovGetEigenval(int m) const { return _covMatrix.GetEigenval(m);} + + //be careful when using this functions!!! + double CovGetEigenval(int m) const { if(_intMatrix.GetEigenval(m)>TOL) return 1./_intMatrix.GetEigenval(m); else return 0.0; } + //const std::vector& CovGetEigenvec(int m) const { return _covMatrix.GetEigenvec(m);} + const std::vector& CovGetEigenvec(int m) const { return _intMatrix.GetEigenvec(m);} + // const std::vector& CovGetEigenvec_fast(int m) const { return _intMatrix.GetEigenvec(m);} + + + void readParameters(const char *filename); + + int readPDBFile(const char* fname); + int readPDBFile(std::string fname){return readPDBFile(fname.c_str());} + int readPDBFile(FILE *fp); + + void setFast(){_FAST=true;} + void setFast(bool fast){_FAST=fast;} + + void setExp(){_EXP=true;} + void setExp(bool fast){_EXP=fast;} + + int getN(){ return _N;} + + void setStructure(Structure3d structure); + void setCutOff(double cutoff){_cutoff=cutoff;} + void setRprot(double r){_R_prot=r;} + void setNtopVectors(int n){_ntop_vectors=n;} + + void setBeadList(vector list){_beadList=list; _structure.setBeadList(_beadList);} + void setOutBeadList(vector list){_outBeadList=list; _TRACCIAMENTO=true;} //ATTENZIONE: qui devo stare attento che non sia outbeadlist>beadlist + + Structure3d & getStructure(){return (_structure);} + + void constructContactMap(); + void constructIntMat(); + void computeCovar(); + int Solve(); + + void Dump(const char *name); + void Dump(std::string name){return Dump(name.c_str());} + void dumpCovEigenvalues(const char*name){_intMatrix.dumpEigenvalues(name);} + void dumpCovEigenvalues(std::string name){return dumpCovEigenvalues(name.c_str());} + // void dumpMSF(const char *name){_intMatrix.dumpMSF(name);} + void dumpMSF(const char *name); + void dumpMSF(std::string name){return dumpMSF(name.c_str());} + void dumpDistFluc(const char *name); + void dumpDistFluc(std::string name){return dumpDistFluc(name.c_str());} + void dumpDistVecFluc(const char *name); + void dumpDistVecFluc(std::string name){return dumpDistVecFluc(name.c_str());} + void dump_top_modes(int, double, const char *name); + void dump_top_modes(int ntop, double amp, std::string name){return dump_top_modes(ntop,amp,name.c_str());} + + double getDistFluc(int, int); + + void dumpNearestNeighbours(const char *name); + void dumpNearestNeighbours(std::string name){return dumpNearestNeighbours(name.c_str());} + + int tracciamento(); //0 if success -1 if not + + void computeForCovMatrix(); + void dumpForCovMatrix(const char*); + void dumpForCovMatrix(std::string name){return dumpForCovMatrix(name.c_str());} + + + private: + // Matrix *_intMatrix; + IntMatrix _intMatrix; + CovMatrix _covMatrix; //ATTENZIONE: posso fare questa cosa? Lo dichiaro senza niente e poi quando lo chiamo gli assegno le cose? + CovMatrix1d _forCovMatrix; + + double _cutoff; + double _R_B; + double _R_BB; + double _R_PP; + double _R_prot; + double _R_N7; + double _Kcov; + double _KBB; + double _KB; + double _KPP; + + int _N; + // int _size; + int _ntop_vectors; + Structure3d _structure ; //, _out_struct; questa volendo la potrei fare. Ma alla fine una volta che traccio mi interessa solo la struttura "ridotta". O no? + + int **_cmap; + + vector _beadList; + + vector _outBeadList; + vector _outResList; + + bool _TRACCIAMENTO; + + bool _COM; + + // bool _PROTEIN; + + bool _OUTRES; // true if you need to reduce the output to a set of residues + + bool _FAST; + + bool _DPF; + + bool _RANDOM; + + bool _EXP; + + bool _RARAND; + + bool _COVFOR; + + bool _DUMPMODES; + bool _DUMPEIGENVALUES; + bool _DUMPCOVMAT; + bool _DUMPDISTFLUC; + bool _DUMPREDCOVMAT; + bool _DUMPNORMREDCOVMAT; + bool _DUMPMSD; + + void constructRandMat(double **); + void constructDemiRandMat(); + void constructExpMat(double **); + + +}; + + + + +#endif /*ELASTICNET_H_*/ diff --git a/EntrComp.C b/EntrComp.C new file mode 100644 index 0000000..1c3d58b --- /dev/null +++ b/EntrComp.C @@ -0,0 +1,138 @@ +/* Program EntrComp: */ +/* reads a pdb trajectory and */ +/* for each frame constructs */ +/* the ENM and compute the entropy */ +/* PhD student at SISSA, Trieste */ +/* March 27th, 2014 */ + + +#include + +#include "Structure3d.h" +#include "ElasticNet.h" +#include + + +#define K_BOLT 1.3806488e-23 // J K^-1 +#define N_AVOG 6.02214129e23 // mol^-1 +#define TEMP 300 // K + +using namespace std; +void help_display(){ + cout<<"Help:"< input file (required) format .pdb"< output file (DEFAULT = entropy.dat)"< parameters file (DEFAULT = PARAMS_RIBOGM.DAT)"< beads to consider (default = C1' C2 P )"< bead_list; + // bead_list.push_back("P"); + // bead_list.push_back("C2"); + // bead_list.push_back("C4'"); + + // ++ READ PARAMETERS + cout<<" --- Program ENTRCOMP --- "<argc-1){cout<<"ERROR: empty parameter"<argc-1){cout<<"ERROR: empty parameter"<argc-1){cout<<"ERROR: empty parameter"< invalid parameter"< +//} + +#define TOL 0.000001 + +#include +#include +#include +#include "omp.h" + + +#include "Matrix.h" +#include "my_malloc.h" +#include "io.h" +//#include "Vector3d.h" + +#include "lapack_matrix_routines_wrapper_v3.cc" + +using namespace std; + +//Constructors for the class Matrix + +Matrix:: Matrix():_DIM(3){ //forse sta cosa non si dovrebbe fare. Se non dichiaro la dimensione della matrice non la creo... + _N=0; + _size=0; + cout<<"Constructing a null matrix"<0){ + _eigenvalues=src._eigenvalues; + _eigenvectors=src._eigenvectors; + } +} + + + +Matrix& Matrix::operator=(Matrix const&src){ //assignment operator + if (this==&src) return *this; + if (_N!=src._N){ + if(_N>0) free_d2t(_matrix); + _N=src._N; + _size=_DIM*_N; + _matrix=d2t(_size,_size); + } + cout<<"Assigning matrix"<0){ + _eigenvalues=src._eigenvalues; + _eigenvectors=src._eigenvectors; + } + return *this; +} + +//############################################# + +void Matrix:: SetN(int N){ + //size of matrix + if(_N==N) return;// AM I SURE THAT THIS IS SAFE?? + if(_size>0) free_d2t(_matrix);// AM I SURE THAT THIS IS SAFE?? + _N=N; + _size=_DIM*_N; + cout<<"Constructing a blank matrix of size N = "<0) free_d2t(_matrix);// AM I SURE THAT THIS IS SAFE?? + _N=N; + _size=_DIM*_N; + cout<<"Matrix:: matrix of size N = "< Matrix::GetEigenvec(int i){ + // vector eigenvec; + // for(int j=0;j<_N;++j){ + // double x=_eigenvectors[i][_DIM*j+0]; + // double y=_eigenvectors[i][_DIM*j+1]; + // double z=_eigenvectors[i][_DIM*j+2]; + // Vector3d ivec(x,y,z); + // eigenvec.push_back(ivec); + // } + // return eigenvec;} + + + +void Matrix::SetEigenvalues(std::vector eval){ + if(_size!=eval.size()){ + cerr<0) _eigenvalues.clear(); + + try{ + _eigenvalues.reserve(eval.size()); + }catch(std::bad_alloc const&){ + cerr << endl << "NormalMode eigenvalues memory allocation failed!" << endl << flush; + exit(1); + } + for(int i=0;i > evec){ + int numModes=evec.size(); + int dim=evec.at(0).size(); + if(_size!=numModes){ + cerr<0) _eigenvectors.clear(); + + try{ + _eigenvectors.resize(numModes); + for (int i = 0; i < numModes; ++i) { + _eigenvectors[i].reserve(dim); + for (int j = 0; j < dim; ++j) { + Vector3d p3d=evec.at(i).at(j); + _eigenvectors[i].push_back(p3d); + } + } + + }catch(std::bad_alloc const&){ + cerr << endl << "NormalMode eigenvector memory allocation failed!" << endl << flush; + exit(1); + } +} + + +void Matrix::Compose(){ + if(_size>0) free_d2t(_matrix);// AM I SURE THAT THIS IS SAFE?? + _size=_eigenvalues.size(); + if(_size<1){ + cout<<"Matrix.Compose --> ERROR: unvalid size of the matrix. Maybe you forgot to set the eigenvalues"< ERROR: unvalid size of the matrix. Maybe you forgot to set the eigenvectors"< eval, std::vector > evec){ + SetEigenvalues(eval); + SetEigenvectors(evec); + _swapEigen();//questo e' giusto solo per una covmatrix!!! + // cout<<"cacca"<0){ + _eigenvalues.clear(); + _eigenvectors.clear(); + cout<<"deleting old eigenstuff"< eigenvec; + // for(int j=0;j<_N;++j){ + // double x=temp_eigenvectors[i][_DIM*j+0]; + // double y=temp_eigenvectors[i][_DIM*j+1]; + // double z=temp_eigenvectors[i][_DIM*j+2]; + // Vector3d ivec(x,y,z); + // eigenvec.push_back(ivec); + // } + // _eigenvalues.push_back(temp_eigenvalues[i]); + // _eigenvectors.push_back(eigenvec); + // } + free_d1t(temp_eigenvalues); + free_d2t(temp_eigenvectors); +} +//########################################### + +void Matrix::dumpAll(const char *name){ + dumpEigenvalues(name); + dumpEigenvectors(name); + dumpFullMatrix(name); + dumpReducedMatrix(name); + dumpNormalizedReducedMatrix(name); +}//END FUNCTION DUMP + + + +void Matrix:: Times(double a){ + /* This function multiplies the matrix for a certain scalar value. + (be carefull. It modifies the matrix itself!!) */ + for(int i=0; i<_size; i++){ + for(int j=0; j<_size; j++){ + _matrix[i][j]*=a; + } + } + if(_eigenvalues.size()>0){ + for(int i=0;i<_size;i++){ + _eigenvalues[i]*=a; + }//enddo i + }//endif +}//end function + +void Matrix:: Plus(const Matrix &M){ + for(int i=0; i<_size; i++){ + for(int j=0; j<_size; j++){ + _matrix[i][j]+=M._matrix[i][j]; + } + } + if(_eigenvalues.size()>0){ + _eigenvalues.clear(); + _eigenvectors.clear(); + cout<<"Deleting eigenstuff: you need to compute them again"<eval_tmp; + vector >evec_tmp; + int nzero=0; + for(int k=0; k <_size; ++k){ + double eval= _eigenvalues[_size-k-1]; + if(fabs(eval) tmp_vec=_eigenvectors[_size-k-1]; + eval_tmp.push_back(tmp_val); + evec_tmp.push_back(tmp_vec); + } + for(int k=0; k tmp_vec=_eigenvectors[_size-k-1]; + eval_tmp.push_back(0.0); + evec_tmp.push_back(tmp_vec); + } //put the zero eigenvalues at the end + tempMat.Compose(eval_tmp,evec_tmp); + cout << "GetInverse: tempo nel ciclo 2 = "<TOL) + logdet+=log(_eigenvalues.at(i)); + } + return logdet; +} + + +//################################################################### +void Matrix::dumpEigenvalues(const char *name){ + char filename[200]; + FILE *fp; + sprintf(filename,"%s_eigenvalues.dat",name); + fp =open_file_w(filename); + for(int i=0; i <_size; i++){ + fprintf(fp,"%4d %e\n",i,_eigenvalues.at(_size-i-1)); + } + fclose(fp); + printf("Eigenvalues written to file %s.\n",filename); +} +/******************************************/ +void Matrix::dumpEigenvectors(const char *name){ + dumpTopVectors(-1,name); +} +/******************************************/ +void Matrix::dumpTopVectors(int ntop,const char *name){ + char filename[200]; + FILE *fp; + if (ntop<0) ntop=_size; + for(int i=0; i < ntop; i++){ + if (i >= _size) break; + sprintf(filename,"%s_eigenvector_%d.dat",name,i); + fp =open_file_w(filename); + for(int j=0; j < _N; j++){ + fprintf(fp,"%4d ",i); + // for(int mu=0; mu < 3; mu++){ + // fprintf(fp,"%lf ",_eigenvectors[_size-i-1][3*j+mu]); + // } + fprintf(fp,"%lf ",_eigenvectors.at(_size-i-1).at(j).X); + fprintf(fp,"%lf ",_eigenvectors.at(_size-i-1).at(j).Y); + fprintf(fp,"%lf ",_eigenvectors.at(_size-i-1).at(j).Z); + + fprintf(fp,"\n"); + } + fclose(fp); + printf("Eigenvector number %4d written to file %s\n",i,filename); + } +} +/******************************************/ +void Matrix::dumpFullMatrix(const char *name){ + char filename[200]; + // FILE *fp; + sprintf(filename,"%s_full_matrix.dat",name); + ofstream fp; + fp.open(filename); + // cout<<"cacca"< eval, std::vector > evec){ + SetEigenvalues(eval); + SetEigenvectors(evec); + // _swapEigen();//questo e' giusto solo per una covmatrix!!! + // cout<<"cacca"< tempvec; + + for(int i=0; i<_size/2; ++i){ + // switch the eigenvalues + tempval=_eigenvalues[i]; + _eigenvalues[i]=_eigenvalues[_size-1-i]; + _eigenvalues[_size-1-i]=tempval; + //switch the eigenvectors + tempvec=_eigenvectors[i]; + _eigenvectors[i]=_eigenvectors[_size-1-i]; + _eigenvectors[_size-1-i]=tempvec; + } +} + +CovMatrix IntMatrix::GetInverse(){ + time_t begin_time2=time(0); + CovMatrix tempMat(_N); + vectoreval_tmp; + vector >evec_tmp; + int nzero=0; + for(int k=0; k <_size; ++k){ + double eval= _eigenvalues[_size-k-1]; + if(fabs(eval) tmp_vec=_eigenvectors[_size-k-1]; + eval_tmp.push_back(tmp_val); + evec_tmp.push_back(tmp_vec); + } + for(int k=0; k tmp_vec=_eigenvectors[_size-k-1]; + eval_tmp.push_back(0.0); + evec_tmp.push_back(tmp_vec); + } //put the zero eigenvalues at the end + // cout<<"%@$%#@$%@#%@#size = "<= _size-6) ivec-=_size; //devo stampare alla fine i 6 vettori nulli! + sprintf(filename,"%s_eigenvector_%d.dat",name,i); + fp =open_file_w(filename); + for(int j=0; j < _N; j++){ + // fprintf(fp,"%4d ",ivec); + // for(int mu=0; mu < 3; mu++){ + // //fprintf(fp,"%lf ",_eigenvectors[i][3*j+mu]); + // fprintf(fp,"%lf ",_eigenvectors[_size-i-1-6][3*j+mu]); + // } + fprintf(fp,"%4d ",ivec); + fprintf(fp,"%lf ",_eigenvectors.at(_size-1-ivec-6).at(j).X); + fprintf(fp,"\n"); + fprintf(fp,"%4d ",ivec); + fprintf(fp,"%lf ",_eigenvectors.at(_size-1-ivec-6).at(j).Y); + fprintf(fp,"\n"); + fprintf(fp,"%4d ",ivec); + fprintf(fp,"%lf ",_eigenvectors.at(_size-1-ivec-6).at(j).Z); + fprintf(fp,"\n"); + } + fclose(fp); + printf("Eigenvector number %4d written to file %s\n",i,filename); + } +} + + +//#################################################################### +//#################################################################### +//#################################################################### + +CovMatrix::CovMatrix():Matrix(){} +CovMatrix::CovMatrix(int i): Matrix(i){} +CovMatrix::CovMatrix(int i, double **d): Matrix(i,d){} + +void CovMatrix::dumpMSF(const char*name){ + char filename[200]; + FILE *fp; + sprintf(filename,"%s_mean_square_displ.dat",name); + fp =open_file_w(filename); + for(int i=0; i < _N; i++){ + double temp=GetMSF(i); + fprintf(fp,"%4d %e\n",i,temp); + } + fclose(fp); + printf("Beads' mean square displacement written to file %s\n",filename); +} + + +double CovMatrix::GetMSF(int i){ + double temp=GetMSF(i,0,0)+GetMSF(i,1,1)+GetMSF(i,2,2); + return temp;} + +double CovMatrix::GetMSF(int i, int mu, int nu){ + double temp= _matrix[_DIM*i+mu][_DIM*i+nu]; + return temp;} + + + + + +IntMatrix CovMatrix::GetInverse(){ + IntMatrix tempMat(_N); + time_t begin_time2=time(0); + vectoreval_tmp; + vector >evec_tmp; + int nzero=0; + for(int k=0; k <_size; ++k){ + double eval= _eigenvalues[k]; + if(fabs(eval) tmp_vec=_eigenvectors[k]; + eval_tmp.push_back(tmp_val); + evec_tmp.push_back(tmp_vec); + } + for(int k=0; k tmp_vec=_eigenvectors[k]; + // eval_tmp.push_back(0.0); + // evec_tmp.push_back(tmp_vec); + eval_tmp.push_back(0.0); + evec_tmp.push_back(tmp_vec); + } //put the zero eigenvalues at the end + tempMat.Compose(eval_tmp,evec_tmp); + cout << "GetInverse: tempo nel ciclo 2 = "< CQ_m_j; + double * CQ_M_j=d1t(_size); + for(int m=0;m<_size;++m){ + for(int l=0;l<_N;++l){ + CQ_M_j[m]+=_matrix[m][_DIM*l+0]*D._eigenvectors[_size-j-1].at(l).X; + CQ_M_j[m]+=_matrix[m][_DIM*l+1]*D._eigenvectors[_size-j-1].at(l).Y; + CQ_M_j[m]+=_matrix[m][_DIM*l+2]*D._eigenvectors[_size-j-1].at(l).Z; + }//enddo l + }//enddo m + for(int i=0;i0){ + _eigenvectors.clear(); + _eigenvalues.clear(); + } + for(int i=0;i<_size;++i){ + for(int j=0;j<_size;++j){ + _matrix[i][j]=temp_mat[i][j]; + } + } + +} + + +double RWSIP(const CovMatrix &A,const CovMatrix &B){ + bool NULL_EVAL=false; + double RWSIP=0.0; + cout<<"### RWSIP ###"< evec_A=A._eigenvectors[m-6]; + for (int n = 6; n < nmodes; ++n) { + double eval_B=B._eigenvalues[n-6]; + if(eval_B evec_B=B._eigenvectors[n-6]; + double scal_prod=0; + for(int i=0; inmodes) return -666; + for (int m = 0; m < ntop; ++m) { + double eval_A=A._eigenvalues[nmodes-1-m]; + if(eval_A evec_A=A._eigenvectors[nmodes-1-m]; + for (int n = 0; n < ntop; ++n){ + double eval_B=B._eigenvalues[nmodes-1-n]; + if(eval_B evec_B=B._eigenvectors[nmodes-1-n]; + double scal_prod=0; + for(int i=0; inmodes) return -666; + for (int m = 0; m < ntop; ++m) { + double eval_A=A._eigenvalues[nmodes-1-m-6]; + if(eval_A evec_A=A._eigenvectors[nmodes-1-m-6]; + for (int n = 0; n < ntop; ++n) { + double eval_B=B._eigenvalues[nmodes-1-n]; + if(eval_B evec_B=B._eigenvectors[nmodes-1-n]; + double scal_prod=0; + for(int i=0; inmodes) return -666; + for (int m = 0; m < ntop; ++m) { + double eval_A=A._eigenvalues[nmodes-1-m-6]; + if(eval_A evec_A=A._eigenvectors[nmodes-1-m-6]; + for (int n = 0; n < ntop; ++n) { + double eval_B=B._eigenvalues[nmodes-1-n-6]; + if(eval_B evec_B=B._eigenvectors[nmodes-1-n-6]; + double scal_prod=0; + for(int i=0; inu_frac*somma) break; + } + cout<<"NTOP="< +#include "Vector3d.h" + +class CovMatrix; +class IntMatrix; + +class Matrix +{ + protected: + int _N; //number of beads; + int _size; //dimension of the matrix + int _DIM; + + double **_matrix; //elements of the correaltion matrix + // TODO Potrebbe servirmi salvare la matrice inversa per non ricalcolarla + //bool _INV; // double **_invmat; + //eigenstuff for the spectral decomposition + std::vector _eigenvalues; + std::vector > _eigenvectors; + + void _swapEigen(); + + public: + Matrix(); + Matrix(int); + Matrix(int, double **); + virtual ~Matrix(); + Matrix(Matrix const &); + + Matrix& operator=(Matrix const &); + + void SetN(int); + void SetMatrix(int, double**); + double& GetElement(int i,int j){return _matrix[i][j];} + + int GetN(){return _N;} + int GetSize(){return _size;} + + double GetEigenval(int i) const {return _eigenvalues.at(i);} + const std::vector& GetEigenvec(int i) const {return _eigenvectors[i];} + + + void SetEigenvalues(std::vector evals); + void SetEigenvectors(std::vector > evec); + + void Compose(); // this function "composes" the matrix of which eigenvalues and eigenvectors are known + void Compose(std::vector evals, std::vector > evec); + + void Decompose(); + + + // void GetInverse(Matrix &outMat); + Matrix GetInverse(); + + double GetTrace(); + double GetLogDet(); + + void Times(double); + void Plus(const Matrix &M); + + void Print(); + void dumpAll(const char *);//better not to use it +//################################################################### + void dumpEigenvalues(const char *name); + void dumpEigenvectors(const char *name); + void dumpTopVectors(int i, const char *name); + + void dumpFullMatrix(const char *name); + void dumpReducedMatrix(const char *name); + void dumpNormalizedReducedMatrix(const char *name); + + void dumpFullInvMatrix(const char *name); + void dumpReducedInvMatrix(const char *name); + + /*********************************/ + //di base stampa sempre le eigenvalues e gli eigenvectors nell'ordine inverso (che per la cov e' da maggiore a minore, per la INT viceversa) + void dumpEigenvalues (std::string name){return dumpEigenvalues (name.c_str());} + void dumpEigenvectors (std::string name){return dumpEigenvectors(name.c_str());} + void dumpTopVectors (int i, std::string name){return dumpTopVectors(i,name.c_str());} + void dumpFullMatrix (std::string name){return dumpFullMatrix (name.c_str());} + void dumpReducedMatrix (std::string name){return dumpReducedMatrix(name.c_str());} + void dumpNormalizedReducedMatrix(std::string name){return dumpNormalizedReducedMatrix (name.c_str());} + void dumpFullInvMatrix (std::string name){return dumpFullInvMatrix (name.c_str());} + void dumpReducedInvMatrix (std::string name){return dumpReducedInvMatrix (name.c_str());} + + //################################################################### + + friend Matrix dot_product(const Matrix &A,const Matrix &B); + friend Matrix dot_product(const IntMatrix &A,const Matrix &B); + friend Matrix dot_product(const Matrix &A,const IntMatrix &B); + friend Matrix dot_product(const IntMatrix &A,const IntMatrix &B); + +}; + + +// NB: nelle CovMatrix gli eigenstuff sono salvati in ordine di eigenvalues crescente! (ma perche' poi?) + + +class IntMatrix : public Matrix { + public: + IntMatrix(); + IntMatrix(int i); + IntMatrix(int i, double **d); + + + friend Matrix dot_product(const IntMatrix &A,const Matrix &B); + friend Matrix dot_product(const Matrix &A,const IntMatrix &B); + friend Matrix dot_product(const IntMatrix &A,const IntMatrix &B); + friend double RWSIP(const IntMatrix &A, const IntMatrix &B); + friend double RWSIP(const IntMatrix &A, const CovMatrix &B); + friend double RWSIP(const CovMatrix &A, const IntMatrix &B){return RWSIP(B,A);} + friend double RMSIP(const IntMatrix &A, const IntMatrix &B, int ntop); + friend double RMSIP(const IntMatrix &A, const CovMatrix &B, int ntop); + friend double RMSIP(const CovMatrix &A, const IntMatrix &B, int ntop){return RMSIP(B,A,ntop);} + + friend double RWSIP_new(const IntMatrix &A, const IntMatrix &B); + friend double RWSIP_new(const IntMatrix &A, const CovMatrix &B); + friend double RWSIP_new(const CovMatrix &A, const IntMatrix &B){return RWSIP_new(B,A);} + + CovMatrix GetInverse(); + void Decompose(); + // void Compose(); + void Compose(std::vector evals, std::vector > evec); + void dumpMSF(const char*); + void dumpMSF(std::string name){return dumpMSF(name.c_str());} + double GetMSF(int i); + double GetMSF(int i, int mu, int nu); + + void dumpCovEigenvalues (std::string name){return dumpCovEigenvalues (name.c_str());} + void dumpCovEigenvectors (std::string name){return dumpCovEigenvectors(name.c_str());} + void dumpTopCovVectors (int i, std::string name){return dumpTopCovVectors(i,name.c_str());} + void dumpCovEigenvalues(const char *name); + void dumpCovEigenvectors(const char *name); + void dumpTopCovVectors(int i, const char *name); + + +}; + + +class CovMatrix : public Matrix { + public: + CovMatrix(); + CovMatrix(int i); + CovMatrix(int i, double **d); + + IntMatrix GetInverse(); + void dumpMSF(const char*); + void dumpMSF(std::string name){return dumpMSF(name.c_str());} + double GetMSF(int i); + double GetMSF(int i, int mu, int nu); + + void Reduce(int NTOP,const CovMatrix &D); + + friend double RWSIP(const CovMatrix &A, const CovMatrix &B); + friend double RWSIP(const IntMatrix &A, const CovMatrix &B); + friend double RWSIP_new(const CovMatrix &A, const CovMatrix &B); + friend double RWSIP_new(const IntMatrix &A, const CovMatrix &B); + friend double RMSIP(const CovMatrix &A, const CovMatrix &B, int ntop); + friend double RMSIP(const IntMatrix &A, const CovMatrix &B, int ntop); + // friend double RWSIP(const CovMatrix &A, const IntMatrix &B){return RWSIP(B,A);} + + friend double ComputeBhatt(CovMatrix A,CovMatrix B); + +}; + + class CovMatrix1d : public CovMatrix { + + protected: + std::vector > _eigenvectors; + void _swapEigen(); + public: + //Costruttori + CovMatrix1d(); + CovMatrix1d(int); + CovMatrix1d(int, double **); + + //func diverse perche' usano gli eigenvectors + const std::vector& GetEigenvec(int i) const {return _eigenvectors[i];} + void SetEigenvectors(std::vector >); + void Compose(); + void Compose(std::vector, std::vector >); + void Decompose(); + Matrix GetInverse(); //TODO + void dumpEigenvectors(const char *); + void dumpTopVectors(int ,const char *); + void dumpMatrix(const char *); + void dumpReducedMatrix(const char *fname) + // {cout<<"*** WARNING: use CovMatrix1d->dumpMatrix() ***"<dumpMatrix() ***"<dumpMatrix() ***"< +#include +#include +#include +#include +#include +#include +#include +#include "NormalModes.h" + +using namespace std; + + + +NormalModes::NormalModes() : _numModes(0){ +} + +NormalModes::~NormalModes() { +} + + +void NormalModes::readFromFile(const std::string& basename, size_t numModes, size_t dim, bool printProgress = false) { + char fname[MAX_PATH_LENGTH] ; + sprintf(fname,"%s%s", basename.c_str(), "_eigenvalues.dat"); + ifstream eigvalFile(fname); + if (!eigvalFile.is_open()){ + cerr << CLRLINE << "Fatal error. Could not open " << fname << endl; + exit(1); + } + + try{ + _eigenvalues.reserve(numModes); + }catch(std::bad_alloc const&){ + cerr << CLRLINE << "NormalMode eigenvalues memory allocation failed!" << endl << flush; + exit(1); + } + + while(eigvalFile.good()){ + if (_eigenvalues.size() == numModes) break; + int index; + double val; + eigvalFile >> index >> val; + _eigenvalues.push_back(val); + } + + if(_eigenvalues.size() < numModes){ + cerr << CLRLINE << "Error. Could read only " << _eigenvalues.size() << " eigenvalues. I wanted "<> index >> p3d; + _eigenvectors[i].push_back(p3d); + + } + if(eigvecFile.fail()){ + cerr << CLRLINE << "Fatal error reading " << fname << endl << flush; + exit(1); + } + if(printProgress) printProgressBar(i,numModes,std::cout); + + } + + }catch(std::bad_alloc const&){ + cerr << CLRLINE << "NormalMode eigenvector memory allocation failed!" << endl << flush; + exit(1); + } + + + _numModes = numModes; + +} + +real_t NormalModes::getEigenvalue(size_t i) const { +assert(i<_eigenvalues.size()); + return _eigenvalues[i]; +} + +const std::vector& NormalModes::getEigenvector(size_t i) const { + assert(i<_eigenvectors.size()); + return _eigenvectors[i]; +} + +void NormalModes::clear() { + _eigenvalues.resize(0); + _eigenvectors.resize(0); + _numModes = 0; +} diff --git a/NormalModes.h b/NormalModes.h new file mode 100644 index 0000000..b51dc68 --- /dev/null +++ b/NormalModes.h @@ -0,0 +1,39 @@ +/* + * NormalModes.h + * + * Created on: Jul 5, 2012 + * Author: gpolles + */ + +#ifndef NORMALMODES_H_ +#define NORMALMODES_H_ +#include +#include "Vector3d.h" +#include "common.h" + + + +class NormalModes +{ +public: + NormalModes(); + virtual ~NormalModes(); + + void readFromFile(const std::string& basename, size_t n, size_t dim, bool printProgress); + + real_t getEigenvalue(size_t i) const; + const std::vector& getEigenvector(size_t i) const; + + size_t getNumModes() const { return _numModes; } + + void clear(); + +private: + std::vector _eigenvalues; + std::vector > _eigenvectors; + size_t _numModes; + +}; + + +#endif /* NORMALMODES_H_ */ diff --git a/PrincipalComp.cc b/PrincipalComp.cc new file mode 100644 index 0000000..5af3844 --- /dev/null +++ b/PrincipalComp.cc @@ -0,0 +1,671 @@ +/* + Created by Giovanni Pinamonti + PhD student @ + Scuola Internazionale Superiore di Studi Avanzati, Trieste, Italy + November 15th 2013 + */ + +//extern "C" { +//#include +//} + + +#include +#include +#include +#include + +#include "PrincipalComp.h" +#include "my_malloc.h" +#include "io.h" +#include "Vector3d.h" + + +using namespace std; + +//######################################## +//Constructors for the class PrincipalComp +//######################################## +PrincipalComp::PrincipalComp():_ntmax(0),_ntmin(0){_N=0;_size=0;} + +PrincipalComp::PrincipalComp(int n):_ntmax(0),_ntmin(0){ + _N=n; _size=3*n;} + +PrincipalComp::PrincipalComp(string fname):_ntmax(0),_ntmin(0){ + readTrajectory(fname);} + +PrincipalComp::PrincipalComp(string fname, std::vector beadlist):_ntmax(0),_ntmin(0){ + setBeadList(beadlist); + readTrajectory(fname.c_str());} + +PrincipalComp::PrincipalComp(string fname, BeadList beadlist):_ntmax(0),_ntmin(0){ + setBeadList(beadlist); + readTrajectory(fname.c_str());} + + +PrincipalComp::~PrincipalComp() { } + +//######################################## + + +void PrincipalComp::readTrajectory(const char* fname){ +// vector list; +// return readTrajectory(fname, list); +// } + +// void PrincipalComp::readTrajectory(const char* fname, std::vector beadlist){ + + Structure3d structure(0); + + structure.setBeadList(_beadList); + + cout<<"PrincipalComp: Opening file"<0)) cout<<"WARNING: old size ("<<_size<<") different from the new one ("<<3*_N<<")"<"<_ntmax)&&(_ntmax>0)) break; + Nsteps_read++; + for(int i=0; i<_N; ++i){ + Vector3d temp_vec_i= structure.getBead(i).getCoordinates(); + for(int mu=0; mu<3; ++mu){ + mean[3*i+mu]+= temp_vec_i[mu]; + } + for(int j=0; j<_N; ++j){ + Vector3d temp_vec_j=structure.getBead(j).getCoordinates(); + for(int mu=0; mu<3; ++mu){ + for(int nu=0; nu<3; ++nu){ + matrix[3*i+mu][3*j+nu]+= temp_vec_i[mu]*temp_vec_j[nu]; + }//enddo nu + }//enddo mu + }//enddo j + }//enddo i + + + //* * * * * * * * * * * * * * * * * * * * * * * * * * * * * + error=structure.readFromPDBFile(fp); // * * * + }//end while // * * * + //* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * + + fclose(fp); + + //normalizing mean and matrix + for(int i=0; i<_size; ++i){ + mean[i]/=Nsteps_read; + for(int j=0; j<_size; ++j){ + matrix[i][j]/=Nsteps_read; + }//enddo + + + }//enddo + //"centering" the covariance matrix + // if (_REF){ + // for(int i=0; i<_size; i++){ + // for(int j=0; j<_size; j++){ + // _matrix[i][j]+=-_mean[i]*_refer[j]-_mean[j]*_refer[i]+_refer[i]*_refer[j]; + // //refer[i]*refer[j]; + // }//enddo + // }//enddo + // } + + for(int i=0; i<_N; ++i){ + double x=mean[3*i+0]; + double y=mean[3*i+1]; + double z=mean[3*i+2]; + Vector3d tempvec(x,y,z); + _mean.getBead(i).setCoordinates(tempvec); + } + + for(int i=0; i<_size; ++i){ + for(int j=0; j<_size; ++j){ + matrix[i][j]-=mean[i]*mean[j]; + }//enddo + }//enddo + + _covMat=(CovMatrix(_N,matrix)); + + free_d1t(mean); + free_d2t(matrix); + + cout<<"...done!"< eigenvalues; + vector > eigenvectors; + cout<<"ATTENZIONE QUI STA LEGGENDO EVAL E EVEC COME SE FOSSERO DA UNA PCA!!!! (giusto?)"<> index >> val; + eigenvalues.push_back(val); + } + if(eigenvalues.size() < numModes){ + cerr << CLRLINE << "Error. Could read only " << eigenvalues.size() << " eigenvalues. I wanted "<> index >> p3d; + eigenvectors[i].push_back(p3d); + } + if(eigvecFile.fail()){ + cerr << CLRLINE << "Fatal error reading " << fname << endl << flush; + exit(1); + } + // if(printProgress) printProgressBar(i,numModes,std::cout); + } + }catch(std::bad_alloc const&){ + cerr << CLRLINE << "NormalMode eigenvector memory allocation failed!" << endl << flush; + exit(1); + } + cout<<"PrincipalComp:ReadFromFile->Compose"< mean_forces; + double *mean_forces=d1t(3*_N); + // for(int i=0;i<_N;i++){ + // mean_forces.push_back(Vector3d(0,0,0)); + // } + vector forces; + while(fgets(line,1000,fp)!=NULL){ + ++nlines; //numero di linee scannate + // cout<<"lines="<>stringa)) continue; //butta in stringa la prima merda di iss, se iss e' vuoto passa alla riga dopo + if (!strncmp(stringa,"natoms=",7)){ //se iss inizia con natoms controlla di avere il giusto numero + iss>>argument; + sscanf(argument,"%d",&nforces); + if((_N!=nforces)&&(_N>0)) cout<<"WARNING: old dimension ("<<_N<<") different from the one of forces file ("<>argument; + + if(!strncmp(argument,"(",1)){ // se dopo f c'e' una parentesi + if(sscanf(argument,"(%dx3)",&temp_nforces)!=1) // allora e' la riga che + cerr<<"WARNING: can't read line"<0){ //se non ho gia' letto forze allora aggiorno la matrice + //$$$ $$$ $$$ $$$ $$$ $$$ $$$ $$$ $$$ $$$ $$$ $$$ $$$ $$$ $$$ + //$$$ $$$ --- Qui aggiorno la matrice (meta') --- $$$ $$$ + if(forces.size()!=nforces) {cerr<<"ERROR: problem with numbers"< _beadList; + BeadList _beadList; + + public: + PrincipalComp(); + PrincipalComp(int n); + PrincipalComp(std::string fname); + PrincipalComp(std::string fname, std::vector beadlist); + PrincipalComp(std::string fname, BeadList beadlist); + virtual ~PrincipalComp(); + + void setNTOP(int n){if(n<0) _NTOP=_size; else _NTOP=n;} + void set_ntmax(int n){_ntmax=n;} + void set_ntmin(int n){_ntmin=n;} + + double CovGetEigenval(int m) const { + return _covMat.GetEigenval(m);} + const std::vector& CovGetEigenvec(int m) const { return _covMat.GetEigenvec(m);} + + CovMatrix & getCovMat(){return (_covMat);} //mica funzionano + IntMatrix & getIntMat(){return (_intMat);} // ste due.... + + void readTrajectory(const char* fname); + void readTrajectory(std::string fname){return readTrajectory(fname.c_str());} + // void readTrajectory(const char* fname, std::vector beadlist); + void readForces(const char*); + void readForces(std::string fname){return readForces(fname.c_str());} + + void Diagonalize(); + + void Dump(const char* fname); + void Dump(string fname){return Dump(fname.c_str());} + void DumpCov(const char* fname); + void DumpCov(string fname){return Dump(fname.c_str());} + + void dumpDistFluc(const char*name); + void dumpMSF(const char*name); + void dumpMSF_partial(const char*name); + + + void dump_top_modes(int,double,const char*); + + void ComputeInt(); + + void DumpInt(const char*fname); + void DumpInt(string fname){return DumpInt(fname.c_str());} + +/* double Converge(double *, int); */ + +/* void SetReference(double *); */ + +//################################################################### +/* void dumpMean(char *name); */ + + void Project(int n_comp, const char *fname){return Project(n_comp,fname,fname);} + void Project(int n_comp, const char *iname, const char *oname ); + + void readFromFile(const char *fname, int nmodes); + void readFromFile(string fname, int nmodes){ + return readFromFile(fname.c_str(),nmodes); } + + void setBeadList(BeadList list){ + _beadList=list;} + void setBeadList(std::vector list){ + _beadList=(BeadList(list));} + +}; diff --git a/RiboGM.C b/RiboGM.C new file mode 100644 index 0000000..e7be970 --- /dev/null +++ b/RiboGM.C @@ -0,0 +1,60 @@ +/* Program RiboGM: */ +/* compute a gaussian network model */ +/* for a given RNA molecule. */ +/* Created by Giovanni Pinamonti */ +/* PhD student at SISSA, Trieste */ +/* November 20th, 2013 */ + + +#include + +#include "Structure3d.h" +#include "ElasticNet.h" + +using namespace std; + + +int main(int argc, char*argv[]){ + + if (argc<1) {cout<<"WARNING: not enough input parameters"< list; + // list.push_back("P"); + // list.push_back("C2"); + // list.push_back("C4'"); + + // structure.setBeadList(list); + // structure.readFromPDBFile(argv[1]); + cout<<"Constructing the gaussian model..."< +#include +#include +#include +#include +#include +#include +//#include "quaternion.h" +//#include "kabsch.cc" +//#include "quaternion.c++" + +/* subroutines per l'allineamento con i quaternioni :) :) :) */ +#include "qcprot.h" +#include "qcprot.cpp" + + +#include "Structure3d.h" +//#include "BeadList.h" + +using namespace std; + +Structure3d::Structure3d() :_size(0),_sizeRes(0), _cutoff(DEFAULT_CUTOFF), _title(""){ +} + + +Structure3d::Structure3d(std::string fname) : _size(0), _sizeRes(0), _cutoff(DEFAULT_CUTOFF), _title(""){ + readFromPDBFile(fname); +} + +Structure3d::Structure3d(double cutoff = DEFAULT_CUTOFF) : _size(0), _sizeRes(0),_cutoff(cutoff), _title(""){ +} +#undef DEFAULT_CUTOFF +Structure3d::Structure3d(double cutoff, std::string fname) : _size(0), _sizeRes(0),_cutoff(cutoff),_title(""){ + readFromPDBFile(fname); +} + +Structure3d::Structure3d(vector beads){ + _beads=beads; + _size=beads.size(); +} + +// Structure3d::Structure3d(int N,double *pos) :_cutoff(cutoff){ +// _size=N +//} + + +Structure3d::~Structure3d() { } + +int Structure3d::readFromPDBFile(const char* fname){ + FILE* fp; + if ((fp = fopen (fname, "r")) == NULL) { + fprintf (stderr,"Could not open file %s.\n.Exiting.\n", fname); + exit (1); } + int idum= readFromPDBFile(fp); + fclose(fp); + return idum; +} + +//return 0 if read succesfully, -1 when file is finished (1 if no atoms in the frame read) +int Structure3d::readFromPDBFile(FILE *fp) { + // For better performance, the vector capacity is preallocated + // and expanded when necessary. + const size_t expandSize = 1000; // size of expansion of capacity + + int discard_alternative_position=0; + int currentChainId = 0; + char line[200], *a; + char field[10]; + double x,y,z,bFactor; + int old_resNum=-999; + + + /* PDB FORMAT + + COLUMNS DATA TYPE FIELD DEFINITION + ------------------------------------------------------------------------------------- + 1 - 6 Record name "ATOM " + + 7 - 11 Integer serial Atom serial number. + + 13 - 16 Atom name Atom name. + + 17 Character altLoc Alternate location indicator. + + 18 - 20 Residue name resName Residue name. + + 22 Character chainID Chain identifier. + + 23 - 26 Integer resSeq Residue sequence number. + + 27 AChar iCode Code for insertion of residues. + + 31 - 38 Real(8.3) x Orthogonal coordinates for X in Angstroms. + + 39 - 46 Real(8.3) y Orthogonal coordinates for Y in Angstroms. + + 47 - 54 Real(8.3) z Orthogonal coordinates for Z in Angstroms. + + 55 - 60 Real(6.2) occupancy Occupancy. + + 61 - 66 Real(6.2) tempFactor Temperature factor. + + 77 - 78 LString(2) element Element symbol, right-justified. + + 79 - 80 LString(2) charge Charge on the atom. + + */ + + + + // ### here I clear the beads to cancel the previous frame (if some) + if(_beads.size()>0) {_beads.clear(); _size=0; _sizeRes=0;} + // TODO: Do I need to reset also the capacity? // + // It will only be usefull if the structure // + // size changes from one frame to the other // + + while (fgets(line, 200, fp)){ + if ((strncmp (line, "ENDMDL", 6)==0)||(strncmp (line, "END", 3)==0)){ + if(_size==0) {cout<<"WARNING: no beads read from file"<0){ + updateNeighbors(); + updateCovBonded();} + return -1; //ritorna -1 se e' finito il file; +}//end function readFromPDBFile() + + +void Structure3d::addBead(Bead temp_bead){ + const size_t expandSize = 1000; // size of expansion of capacity + if(_beadList.check(temp_bead)){ + // bool cacca=_beadList.check(temp_bead); + // cout<_size-1) continue; + if (_beads[j].getAtomType()=="P"){ + double dsq=_beads[j].getCoordinates().distanceSQ(_beads[i].getCoordinates()); + if (dsq < CHAIN_BREAK_DIST_SQ){ + _beads[i].addCovBond(&_beads[j]); + _beads[j].addCovBond(&_beads[i]); + } + }//endif + }//endif C4' + }//enddo I +} + +bool Structure3d::isBond(int i, int j){ + for(int k=0;k<_beads[i].getNumCovBonded();++k){ + int n=(*(_beads[i].getCovBonded().at(k))).getIndex(); + if(n==j) return true; + } + return false; +} + +// bool Structure3d::_isABead(char *a){ + +// bool cacca=false; + +// if(_beadList.size()==0) return true; + +// stringstream ss; +// string s; +// ss<>s; + +// for(int i=0;i<_beadList.size();++i) +// { +// string b; +// b=_beadList.at(i); +// if (s.compare(b)==0) cacca=true; +// // cout< sugar, base; + + + // ### here I clear the beads to cancel the previous frame (if some) + if(_beads.size()>0) _beads.clear(); + + while (fgets(line, 200, fp)){ + if (strncmp (line, "ENDMDL", 6)==0){ + updateNeighbors(); + updateCovBonded(); + return 0; + } + /* check if line buffer begins with "ATOM" */ + if ( (!strncmp (line, "ATOM", 4))||(!strncmp (line, "HETATM", 6)) ){ + a = line +12; /* advance line buffer and check the atom's type */ + /* check if char at 16 is diff from ' ' or 'A' */ + discard_alternative_position=1; + if (line[16]==' ') discard_alternative_position=0; + if (line[16]=='A') discard_alternative_position=0; + + if (discard_alternative_position==1) continue; + + a = line +12; /* advance line buffer and check if + we have a suitable atom */ + + bool civa=false; + int resNum; + stringstream ss; + string s, type; + ss<>s; + + /* Phosphates */ + if (s.compare("P")==0){ civa=true; + type=s; + /* Residue number */ + a=line+22; + strncpy(field,a,4); field[4]='\0'; + if (sscanf(field,"%d",&resNum)!=1){fprintf(stderr,"Error. Could not read residue number in molecule.\nLine %s\n",line); exit(1);} + /* Coordinates */ + a=line+30; + strncpy(field,a,8);field[8]='\0'; + if (sscanf(field,"%lf",&x)!=1){fprintf(stderr,"Fatal error. Could not read bead coordinate in molecule.\nLine %s\n",line); exit(1);} + a=line+38; + strncpy(field,a,8);field[8]='\0'; + if (sscanf(field,"%lf",&y)!=1){fprintf(stderr,"Fatal error. Could not read bead coordinate in molecule.\nLine %s\n",line); exit(1);} + a=line+46; + strncpy(field,a,8);field[8]='\0'; + if (sscanf(field,"%lf",&z)!=1){fprintf(stderr,"Fatal error. Could not read bead coordinate in molecule.\nLine %s\n",line); exit(1);} + /* B factor */ + bFactor=0.0; + } + + /* sugars */ + bool issugar=false; + if(s.compare("C2'")==0) issugar=true; + if(s.compare("C3'")==0) issugar=true; + if(s.compare("C4'")==0) issugar=true; + if(s.compare("O4'")==0) issugar=true; + if(s.compare("C1'")==0) issugar=true; + + if(issugar){ + /* store coordinates of sugar atoms */ + double tempx, tempy, tempz; + + a=line+30; strncpy(field,a,8);field[8]='\0'; + if (sscanf(field,"%lf",&tempx)!=1){fprintf(stderr,"Fatal error. Could not read bead coordinate in molecule.\nLine %s\n",line); exit(1);} + + a=line+38; strncpy(field,a,8);field[8]='\0'; + if (sscanf(field,"%lf",&tempy)!=1){fprintf(stderr,"Fatal error. Could not read bead coordinate in molecule.\nLine %s\n",line); exit(1);} + + a=line+46; strncpy(field,a,8);field[8]='\0'; + if (sscanf(field,"%lf",&tempz)!=1){fprintf(stderr,"Fatal error. Could not read bead coordinate in molecule.\nLine %s\n",line); exit(1);} + + sugar.push_back(Vector3d(tempx,tempy,tempz)); + + /* if the ring is completely stored, compute the center */ + if(sugar.size()==5) + { civa=true; bFactor=0.0; type="S"; + x=0;y=0;z=0; + for (int i=0; i<5;++i){ + x+=sugar.at(i).X; + y+=sugar.at(i).Y; + z+=sugar.at(i).Z; + } + x*=0.2; + y*=0.2; + z*=0.2; + sugar.clear(); + /* Residue number */ a=line+22; strncpy(field,a,4); field[4]='\0'; + if (sscanf(field,"%d",&resNum)!=1){fprintf(stderr,"Error. Could not read residue number in molecule.\nLine %s\n",line); exit(1);} } + + }//end issugar + + + /* bases */ + + bool isbase=false; + if(s.compare("C5")==0) isbase=true;// x=0; y=0; z=0;} + if(s.compare("C2")==0) isbase=true; + + if(isbase){ + /* store coordinates of base atoms */ + double tempx, tempy, tempz; + + a=line+30; strncpy(field,a,8);field[8]='\0'; + if (sscanf(field,"%lf",&tempx)!=1){fprintf(stderr,"Fatal error. Could not read bead coordinate in molecule.\nLine %s\n",line); exit(1);} + + a=line+38; strncpy(field,a,8);field[8]='\0'; + if (sscanf(field,"%lf",&tempy)!=1){fprintf(stderr,"Fatal error. Could not read bead coordinate in molecule.\nLine %s\n",line); exit(1);} + + a=line+46; strncpy(field,a,8);field[8]='\0'; + if (sscanf(field,"%lf",&tempz)!=1){fprintf(stderr,"Fatal error. Could not read bead coordinate in molecule.\nLine %s\n",line); exit(1);} + + base.push_back(Vector3d(tempx,tempy,tempz)); + + /* if the ring is completely stored, compute the center */ + if(base.size()==2){ + civa=true; bFactor=0.0; type="B"; + x=0;y=0;z=0; + for (int i=0; i<2;++i){ + x+=base.at(i).X; + y+=base.at(i).Y; + z+=base.at(i).Z; + } + x*=0.5; + y*=0.5; + z*=0.5; + base.clear(); + /* Residue number */ a=line+22; strncpy(field,a,4); field[4]='\0'; + if (sscanf(field,"%d",&resNum)!=1){fprintf(stderr,"Error. Could not read residue number in molecule.\nLine %s\n",line); exit(1);} + } + }//end isbase + + + /**********************************/ + if (civa) { + /* check capacity */ + if(_beads.capacity() <= _beads.size()) _beads.reserve(_beads.capacity() + expandSize); + /* add record */ + _beads.push_back(Bead()); + _beads.back().setCoordinates(Vector3d(x,y,z)); x=0; y=0; z=0; + _beads.back().setBetaFactor(bFactor); + _beads.back().setPdbPosition(_size); + _beads.back().setIndex(_size); + _beads.back().setAtomType(type); + _beads.back().setResNum(resNum); + ++_size; + if (type=="B") ++_sizeRes; + }//endif (civa) + }//endif ATOM + }//END WHILE + if(_beads.size()>0){ + updateNeighbors(); + updateCovBonded();} + return -1; +} + + + + + int Structure3d::dumpPDBFile(const char* fname){ + char name[200]; + ofstream dumpfile; + sprintf(name,"%s",fname); + dumpfile.open(name); + return dumpPDBFile(dumpfile); + } + +int Structure3d::dumpPDBFile(std::ofstream &dumpfile){ + char line[200], *a; + char field[10]; + + + // PDB FORMAT + +// COLUMNS DATA TYPE FIELD DEFINITION + + char record[6]; // 1 2 3 4 5 6 Record name "ATOM " + char serial[5]; // 7 8 9 10 11 Integer serial Atom serial number. + char name[4]; // 13 14 15 16 Atom name Atom name. + char altLoc[2]; // 17 Character altLoc Alternate location. + char resName[3]; // 18 19 20 Residue name resName Residue name. + char chainID[2]; // 22 Character chainID Chain identifier. + char resSeq[4]; // 23 24 25 26 Integer resSeq Residue sequence. + char iCode[2]; // 27 AChar iCode Code 4 insertion of res + char X[7]; // 31 32 33 34 35 36 37 38 Real(8.3) x X in Angstroms + char Y[7]; // 39 40 41 42 43 44 45 46 Real(8.3) y Y in Angstroms + char Z[7]; // 47 48 49 50 51 52 53 54 Real(8.3) z Z in Angstroms + char occupancy[6]; // 55 56 57 58 59 60 Real(6.2) occupancy Occupancy. + char tempFactor[6]; // 61 62 63 64 65 66 Real(6.2) tempFactor B-factor. + char element[2]; // 77 78 LString(2) element Element symbol, right-justified. + char charge[2]; // 79 80 LString(2) charge Charge on the atom. + + + dumpfile<<"REMARK GENERATED FROM STRUCTURE3D"< the two structures to align doesn't have the same size!"< + +#include "Vector3d.h" +#include "Bead.h" +#include +#include "my_malloc.h" +#include "BeadList.h" + + +class Structure3d +{ + + public: + Structure3d(); + Structure3d(double cutoff); + Structure3d(double cutoff, std::string fname); + Structure3d(std::string fname); + + Structure3d(std::vector); + /* Structure3d(int,double*); */ + virtual ~Structure3d(); + + + std::string getTitle(){ return _title;} + void setTitle(std::string title){_title=title;} + void setTitle(const char* title){ std::string str_title(title); return setTitle(str_title);} + + int readFromPDBFile(FILE* fp); + // int readFromPDBFile(std::ifstream $fp); + int readFromPDBFile(const char* fname); + int readFromPDBFile(std::string fname) {return readFromPDBFile(fname.c_str());} + + int centersFromPDBFile(FILE* fp); + int centersFromPDBFile(const char* fname); + int centersFromPDBFile(std::string fname) {return centersFromPDBFile(fname.c_str());} + + int dumpPDBFile(std::ofstream &dumpfile); + int dumpPDBFile(const char* fname); + int dumpPDBFile(std::string fname) {return dumpPDBFile(fname.c_str());} + + void addBead(Bead ); + + + double getDistance(size_t i, size_t j) const; + double getDistanceSQ(size_t i, size_t j) const; + Vector3d getDistanceVector(size_t i, size_t j); + size_t getSize() const {return _size;} + size_t getSizeRes() const {return _sizeRes;} + void updateNeighbors(); + void updateCovBonded(); + + std::vector& getBeads() { + return _beads; + } + Bead& getBead(size_t i){ + return _beads[i]; + } + + + void setBeadList(vector list) + { + _beadList=(BeadList(list)); + } + + void setBeadList(BeadList list) + { + _beadList=list; + } + + + bool isBond(int, int); + + + + double fit(Structure3d *ref_struc, double *w); + + private: + std::vector _beads; + std::vector _coordinates; + std::vector _bFactors; + + std::string _title; + + size_t _sizeRes; + size_t _size; + double _cutoff; + + // bool _isABead(char *); + + // std::vector _beadList; + BeadList _beadList; + +}; + +#endif /* STRUCTURE3D_H_ */ diff --git a/Vector3d.h b/Vector3d.h new file mode 100644 index 0000000..04292aa --- /dev/null +++ b/Vector3d.h @@ -0,0 +1,129 @@ +/* + * Vector3d.h + * + * Created on: Jul 2, 2012 + * Author: gpolles + */ + +#ifndef POINT3D_H_ +#define POINT3D_H_ + +#define SQ(i) (i)*(i) +#include +#include +#include + +class Vector3d +{ + friend std::ostream& operator <<(std::ostream &out, const Vector3d& p); + friend std::istream& operator >>(std::istream &in, Vector3d& p); + + +public: + double X; + double Y; + double Z; + public: + Vector3d(){X=0;Y=0;Z=0;} + Vector3d(double ax, double ay, double az) {X=ax;Y=ay;Z=az;} + Vector3d(const Vector3d& a) {X = a.X ; Y = a.Y ; Z = a.Z;} + + double & getCoord(int mu){ + if(mu==0) return X; + if(mu==1) return Y; + if(mu==2) return Z; + } + + + + bool operator ==(const Vector3d& a) { if(a.X == X && a.Y==Y && a.Z==Z) return true; return false;} + inline Vector3d operator-(const Vector3d& a) const{ + return Vector3d(X-a.X,Y-a.Y,Z-a.Z); + } + inline Vector3d operator-() const{ + return Vector3d(-X,-Y,-Z); + } + inline Vector3d operator+(const Vector3d& a) const{ + return Vector3d(X+a.X,Y+a.Y,Z+a.Z); + } + inline Vector3d operator/(const Vector3d& a) const{ + return Vector3d(X/a.X,Y/a.Y,Z/a.Z); + } + inline Vector3d& operator-=(const Vector3d& a) { + X-=a.X; Y-=a.Y; Z-=a.Z; + return *this; + } + inline Vector3d& operator/=(const double a) { + X/=a; Y/=a; Z/=a; + return *this; + } + inline Vector3d& operator*=(const double a) { + X*=a; Y*=a; Z*=a; + return *this; + } + inline Vector3d& operator+=(const Vector3d& a) { + X+=a.X; Y+=a.Y; Z+=a.Z; + return *this; + } + Vector3d operator*(const double a) const{ + return Vector3d(a*X,a*Y,a*Z); + } + inline Vector3d operator/(const double a) const{ + return Vector3d(X/a,Y/a,Z/a); + } + + double operator[](size_t i) { + return *(&X+i); + } + double operator[](int i) { + return *(&X+i); + } + + + //prodotti scalare e vettoriale (ho cambiato i nomi da dotProducts a dot) + inline double dot(const Vector3d& p) const {return X*p.X + Y*p.Y + Z*p.Z;} + inline Vector3d cross(const Vector3d& p) const { + return Vector3d( Y*p.Z - Z*p.Y, + Z*p.X - X*p.Z, + X*p.Y - Y*p.X); + } + + + inline double distance(const Vector3d& p) const {return sqrt(SQ(p.X-X)+SQ(p.Y-Y)+SQ(p.Z-Z)); } //distanza fra due punti + inline double distanceSQ(const Vector3d& p) const {return SQ(p.X-X)+SQ(p.Y-Y)+SQ(p.Z-Z); } //distanza quadra + + inline double normSQ() const {return X*X + Y*Y + Z*Z;} //norma/modulo quadra/o + inline double norm() const {return sqrt(X*X + Y*Y + Z*Z);} //norma/modulo + + inline double angle(const Vector3d& p) const { + double x = dot(p)/norm()/p.norm(); + return acos( x );} +}; + +//######################################################## +//### funzioni utili che utilizzano la suddetta classe ### +//######################################################## + +inline double norm(Vector3d& v) {return sqrt(v.X*v.X + v.Y*v.Y + v.Z*v.Z);} +inline double normSQ(Vector3d& v) {return (v.X*v.X + v.Y*v.Y + v.Z*v.Z);} + +inline std::ostream& operator <<(std::ostream &out, const Vector3d& p){ + out << "(" << p.X << ", " << p.Y << ", " << p.Z << ")"; + return out; +} + +inline std::istream& operator >>(std::istream &in, Vector3d& p){ + in >> p.X >> p.Y >> p.Z; + return in; +} + +inline Vector3d operator *(const double x, const Vector3d& v) { + return v*x; +} + + +inline double dot(const Vector3d& p,const Vector3d& q) { + return q.X*p.X + q.Y*p.Y + q.Z*p.Z; +} + +#endif /* POINT3D_H_ */ diff --git a/ang.C b/ang.C new file mode 100644 index 0000000..d616265 --- /dev/null +++ b/ang.C @@ -0,0 +1,499 @@ +/* written by Giovanni Pinamonti, October 2014 */ +//============================================================================ +// originally based on +// Name : generateConfigs.cpp +// Author : Guido Polles +// Version : +// Copyright : +// Description : +//============================================================================ + +#include +#include +#include +#include +#include "omp.h" +#include + +#include "Structure3d.h" +#include "PrincipalComp.h" +#include "ElasticNet.h" +//#include "Matrix.h" +using namespace std; + +void help_display(){ + cout<<"Help:"< input file (required) format .pdb"< normalized temperature (default =0.1)"< number of frames (default = -1)"< number of modes to consider (default = all)"< parameters file (DEFAULT = PARAMS_RIBOGM.DAT)"< PCA file (DEFAULT = none) format .pdb"< output file for fluctuations (DEFAULT = none)"< output file for correlations (DEFAULT = none)"< input file with shape data [data in third column] (DEFAULT = none)"< compute the Bhattacharya distance (DEFAULT = false)"< is a riboswitch? (do I need to take away the last res because it's an ADA?) (default=false)"< dump the ENM results on file (default=false)"< are there pre-constructed ENM files? (DEFAULT = false)"< (DEFAULT = intput file name)"< bead_list; + bool out_has_name=false; + bool shape_has_name=false; + bool corr_has_name=false; + bool pca_has_name=false; + bool comp_bhatt=false; + bool ribo=false; + bool dump=false; + bool make_enm=true; + bool has_name=false; + // bead_list.push_back("C1'"); + // bead_list.push_back("C2"); + // bead_list.push_back("P"); + + cout<<" --- Program GenConf --- "<argc-1){cout<<"ERROR: empty parameter"<argc-1){cout<<"ERROR: empty parameter"<argc-1){cout<<"ERROR: empty parameter"<argc-1){cout<<"ERROR: empty parameter"<argc-1){cout<<"ERROR: empty parameter"<argc-1){cout<<"ERROR: empty parameter"<argc-1){cout<<"ERROR: empty parameter"<argc-1){cout<<"ERROR: empty parameter"<argc-1){cout<<"ERROR: empty parameter"<argc-1){cout<<"ERROR: empty parameter"< invalid parameter"< nd(0.0,1); // definisco la distribuzione normale da cui estrarro i numeri gaussiani + boost::variate_generator > rnd(rng, nd); + + if((NTOP<1)||(NTOP>nmodes-6)) NTOP=nmodes-6; + vector sqrt_eval; + for (int m = 6; m < NTOP+6; ++m) { + double s=ENM.CovGetEigenval(nmodes-m-1); + if(s<0.0){cout<<"genConf *** WARNING: negative eigenvalue -> "< coords(n_beads); + + int NFPRINT=NFRAMES/10; + if(NFPRINT<1) NFPRINT=1; + for (size_t f = 0; f < NFRAMES; ++f){ + if(f%NFPRINT==0) cout< width; + for (int m = 6; m < NTOP+6; ++m) { + //double s=ENM.CovGetEigenval(nmodes-m-1); + double s=sqrt_eval[m-6]; + width.push_back(rnd()*s); + } + // $ $ $ $ $ $ $ $ +#pragma omp parallel for + for (size_t i = 0; i < n_beads; ++i){ + for (int m = 6; m < NTOP+6; ++m) { + coords[i] += ENM.CovGetEigenvec(nmodes-m-1)[i]*width[m-6]; + } + } + // $ $ $ $ $ $ $ $ + // # # # # # # # + for (size_t i = 0; i < indexC2.size(); ++i){ //### qui calcolo le distanze + double temp_dsq=coords[i].distanceSQ(coords[i+1]); + d_sq_mean[i]+=temp_dsq; + dist_mean[i]+=sqrt(temp_dsq); + } + + }//enddo frames + + for (size_t i = 0; i < n_C2; ++i){ + // for (size_t i = 0; i < indexC2.size()-2; ++i){ //NB: -2 per escludere ADA del riboswitch, altrimenti -1 !! + d_sq_mean[i]/=NFRAMES; + dist_mean[i]/=NFRAMES; + fluct.push_back(d_sq_mean[i]-(dist_mean[i]*dist_mean[i])); + } + }//endelse + + /// +++ OPEN the OUTPUT file +++ + if(out_has_name){ + ofstream fout; + fout.open(oname); + // ENM.setFast(false); + // ENM.computeCovar(); + + // ENM.Dump("prova"); + // ENM.dump_top_modes(10,5,"prova"); + // ### qui calcolo le distanze e le printo ### + // ENM.dumpDistFluc("prova"); + for (size_t i = 0; i < n_C2; ++i){ + // for (size_t i = 0; i < indexC2.size()-2; ++i){ //NB: -2 per escludere ADA del riboswitch, altrimenti -1 !! + fout< shape; + ifstream shapefile(shape_name); + char cdum[2]; + int idum; + double ddum; + double tmpshape; + /// while(shapefile >> cdum >> idum >> tmpshape >> cdum >> idum >>ddum >> cdum >> idum >>ddum >> cdum >> idum >>ddum >> cdum >> idum >>ddum){ + string line; + while(getline(shapefile,line)){ + std::istringstream iss(line); + iss >> cdum >> idum >> tmpshape; + if(!strncmp(cdum,"#",1)) continue; + shape.push_back(tmpshape); + } + cout<<"### corr ###"<0.0) + C+=1; + else + if((fluct[i]-fluct[j])*(shape[i]-shape[j])<0.0) + D+=1; + }//enddo j + }//enddo i + kend=(C-D)/C_n_2; + corr=(sum_XY-(sum_X*sum_Y)/ndata)/(sqrt((sum_X2-(sum_X*sum_X)/ndata)*(sum_Y2-(sum_Y*sum_Y)/ndata))); + cout< evec_ENM=ENM.CovGetEigenvec(m-6); + for (int n = 6; n < nmodes; ++n) { + double eval_PCA=PCA.getCovMat().GetEigenval(n); + vector evec_PCA=PCA.getCovMat().GetEigenvec(n); + double scal_prod=0; + for(int i=0; inu_frac*somma) break; + } + cout<<"NTOP="< +#include +#include +#include +#include "omp.h" +#include + +#include "Structure3d.h" +#include "PrincipalComp.h" +#include "ElasticNet.h" +//#include "Matrix.h" +using namespace std; + +void help_display(){ + cout<<"Help:"< input file (required) format .pdb"< normalized temperature (default =0.1)"< number of frames (default = -1)"< number of modes to consider (default = all)"< parameters file (DEFAULT = PARAMS_RIBOGM.DAT)"< PCA file (DEFAULT = none) format .pdb"< output file for fluctuations (DEFAULT = none)"< output file for correlations (DEFAULT = none)"< input file with shape data [data in third column] (DEFAULT = none)"< is a riboswitch? (do I need to take away the last res because it's an ADA?) (default=false)"< dump the ENM results on file (default=false)"< are there pre-constructed ENM files? (DEFAULT = false)"< (DEFAULT = intput file name)"< bead_list; + bool out_has_name=false; + bool shape_has_name=false; + bool corr_has_name=false; + bool pca_has_name=false; + bool comp_bhatt=false; + bool ribo=false; + bool dump=false; + bool make_enm=true; + bool has_name=false; + // bead_list.push_back("C1'"); + // bead_list.push_back("C2"); + // bead_list.push_back("P"); + + cout<<" --- Program ang_C1p_P --- "<argc-1){cout<<"ERROR: empty parameter"<argc-1){cout<<"ERROR: empty parameter"<argc-1){cout<<"ERROR: empty parameter"<argc-1){cout<<"ERROR: empty parameter"<argc-1){cout<<"ERROR: empty parameter"<argc-1){cout<<"ERROR: empty parameter"<argc-1){cout<<"ERROR: empty parameter"<argc-1){cout<<"ERROR: empty parameter"<argc-1){cout<<"ERROR: empty parameter"<argc-1){cout<<"ERROR: empty parameter"< invalid parameter"< index_C1p_P; + for (size_t i = 0; i < ref_struc.getSize(); ++i){ + if(ref_struc.getBead(i).getHet()) continue; + if ((ref_struc.getBead(i).getAtomType()=="C1'")||(ref_struc.getBead(i).getAtomType()=="P")||(ref_struc.getBead(i).getAtomType()=="C2")) {is_C1p_P[i]=1; index_C1p_P.push_back(i);} + } + cout<<"N. beads = "< fluct; //vettore su cui salvero' le fluttuazioni C2-C2 + + if(NFRAMES<3){ cout<<"ERRORRE"< nd(0.0,1); // definisco la distribuzione normale da cui estrarro i numeri gaussiani + boost::variate_generator > rnd(rng, nd); + + if((NTOP<1)||(NTOP>nmodes-6)) NTOP=nmodes-6; + vector sqrt_eval; + for (int m = 6; m < NTOP+6; ++m) { + double s=ENM.CovGetEigenval(nmodes-m-1); + if(s<0.0){cout<<"genConf *** WARNING: negative eigenvalue -> "< coords(index_C1p_P.size()); + + // +++ OPEN the OUTPUT file +++ + char oname[1024]; + ofstream dumpfile; + sprintf(oname,"%s_trajectory.pdb",file_name); + dumpfile.open(oname); + // +++ +++ +++ +++ +++ +++ +++ + + int NFPRINT=NFRAMES/10; + if(NFPRINT<1) NFPRINT=1; + for (size_t f = 0; f < NFRAMES; ++f){ + if(f%NFPRINT==0) cout< width; + for (int m = 6; m < NTOP+6; ++m) { + //double s=ENM.CovGetEigenval(nmodes-m-1); + double s=sqrt_eval[m-6]; + width.push_back(rnd()*s); + } + // $ $ $ $ $ $ $ $ +#pragma omp parallel for + for (size_t i = 0; i < index_C1p_P.size(); ++i){ + int i_bead=index_C1p_P.at(i); + for (int m = 6; m < NTOP+6; ++m) { + coords[i] += ENM.CovGetEigenvec(nmodes-m-1)[i_bead]*width[m-6]; //sbagliavo qui... + } + new_struc.getBead(i).setCoordinates(coords[i]); + } + new_struc.dumpPDBFile(dumpfile); + // $ $ $ $ $ $ $ $ + + // # # # # # # # + for (size_t i = 0; i < index_C1p_P.size()-3; ++i){ //### qui calcolo le distanze + if (ref_struc.getBead(index_C1p_P.at(i)).getAtomType()=="C1'"){ + Vector3d d12=coords[i]-coords[i+2]; + Vector3d d23=coords[i+3]-coords[i+2]; + double coseno=d12.dot(d23); + coseno/=d12.norm(); + coseno/=d23.norm(); + double angle=acos(coseno); + an_sq_mean[i]+=angle*angle; + angle_mean[i]+=angle; + } + } + + }//enddo frames + + dumpfile.close(); + for (size_t i = 0; i < n_C1p_P-3; ++i){ + an_sq_mean[i]/=NFRAMES; + angle_mean[i]/=NFRAMES; + // if (ref_struc.getBead(i).getAtomType()=="C1'"){ + fluct.push_back(an_sq_mean[i]-(angle_mean[i]*angle_mean[i])); + // } + } + }//endelse + + /// +++ OPEN the OUTPUT file +++ + if(out_has_name){ + ofstream fout; + fout.open(oname); + // int idum=0; + for (size_t i = 0; i < n_C1p_P-3; ++i){ + if(ref_struc.getBead(index_C1p_P.at(i)).getAtomType()=="C1'"){ + fout< shape; + ifstream shapefile(shape_name); + char cdum[2]; + int idum; + double ddum; + double tmpshape; + string line; + while(getline(shapefile,line)){ + std::istringstream iss(line); + iss >> cdum >> idum >> tmpshape; + if(!strncmp(cdum,"#",1)) continue; + shape.push_back(tmpshape); + } + cout<<"### corr ###"<0.0) + C+=1; + else + if((fluct[i]-fluct[j])*(shape[i]-shape[j])<0.0) + D+=1; + }//enddo j + }//enddo i + kend=(C-D)/C_n_2; + corr=(sum_XY-(sum_X*sum_Y)/ndata)/(sqrt((sum_X2-(sum_X*sum_X)/ndata)*(sum_Y2-(sum_Y*sum_Y)/ndata))); + cout< +#include +#include +#include +#include + + +typedef double real_t; +typedef int int_t; +typedef unsigned int uint_t; + +#define CLRLINE "\r \r" + + +inline void printProgressBar(long double stat, long double max, std::ostream& out, int lenght = 40){ + if (stat>max) stat = max; + int nch = static_cast(floor((stat/max)*lenght)); + out.setf(std::ios::fixed); + out << "\r[" ; + for (int i = 0; i < 40; ++i) { + if (nch < i) out << " "; + else out << "#"; + } + out << "] " << std::setw(6) << std::setprecision(2) + << (stat/max)*100.0 << "%" << std::flush; +} + + +inline std::string str(double x) { + std::stringstream ss; + ss << x; + return ss.str(); +} + +#endif /* COMMON_H_ */ diff --git a/compdist.C b/compdist.C new file mode 100644 index 0000000..6420e2d --- /dev/null +++ b/compdist.C @@ -0,0 +1,167 @@ +/* Program written by Giovanni Pinamonti + PhD student at + Scuola Internazionale Superiori di Studi Avanzati, Trieste, Italy + begin june 11th 2013 */ + +#include +#include +#include +#include +#include +#include +#include + +#include "Matrix.h" //definition of the class CMatrix +#include "Structure3d.h" + +using namespace std; + + +void help_display(){ + cout<<"Help:"< input trajectory file (required) format .pdb"< number of frames (default = 10000)"< beads to consider (default = [ C1' C2 P ] )"< output file for fluctuations (DEFAULT = none)"< is a riboswitch? (do I need to take away the last res because it's an ADA?) (default=false)"<argc-1){cout<<"ERROR: empty parameter"<argc-1){cout<<"ERROR: empty parameter"< invalid parameter"< indexC2; + for (size_t i = 0; i < ref_struc.getSize(); ++i){ + if(ref_struc.getBead(i).getHet()) continue; + if ((ref_struc.getBead(i).getAtomType()=="C2") || + (ref_struc.getBead(i).getAtomType()=="C6")) + {is_C2[i]=1; indexC2.push_back(i);}// cout<<"eccolo: "< coords(indexC2.size()); + + dist_mean=d1t(n_beads); + d_sq_mean=d1t(n_beads); + + cout<<"PrincipalComp: reading trajectory..."< fluct; + int n_C2=0; + if(ribo) n_C2=indexC2.size()-2;//NB: -2 per escludere ADA del riboswitch, altrimenti -1 !! + else n_C2=indexC2.size()-1; + + for (size_t i = 0; i < n_C2; ++i){ + d_sq_mean[i]/=Nsteps; + dist_mean[i]/=Nsteps; + fluct.push_back(d_sq_mean[i]-(dist_mean[i]*dist_mean[i])); + } + /// +++ OPEN the OUTPUT file +++ + // if(out_has_name){ + ofstream fout; + fout.open(oname); + // ### qui calcolo le distanze e le printo ### + for (size_t i = 0; i < n_C2; ++i){ //NB: -2 per escludere ADA del riboswitch, altrimenti -1 !! + fout<<"X "< + +double correlation(std::vector data1, std::vector data2){ + double corr=0.0; + int ndata1=data1.size(); + if(!(ndata1==data2.size())){std::cerr<<"ERROR: correlation->lunghezza dati diverse. Non posso correlarle. "< data1, std::vector data2){ + int ndata=data1.size(); + double kend=0.0; + if(!(ndata==data2.size())){cout<<"ERROR: lunghezza dati data1 diversa da fluttuazioni. Non posso correlarle. "<0.0) + C+=1; + else + if((data2[i]-data2[j])*(data1[i]-data1[j])<0.0) + D+=1; + }//enddo j + }//enddo i + kend=(C-D)/C_n_2; + } + return kend; +} + diff --git a/countn.C b/countn.C new file mode 100644 index 0000000..2b49d4b --- /dev/null +++ b/countn.C @@ -0,0 +1,87 @@ +#include +#include +#include +#include +#include "omp.h" +#include + + + +#include "Structure3d.h" +#include "PrincipalComp.h" +#include "ElasticNet.h" +//#include "Matrix.h" +using namespace std; + + +void help_display(){ + cout<<"Help:"< input file (required) format .pdb"< parameters file (DEFAULT = PARAMS_RIBOGM.DAT)"< output name"<argc-1){cout<<"ERROR: empty parameter"<argc-1){cout<<"ERROR: empty parameter"<argc-1){cout<<"ERROR: empty parameter"< invalid parameter"< +#include +#include +#include +#include "omp.h" +#include + +#include "Structure3d.h" +#include "PrincipalComp.h" +#include "ElasticNet.h" +//#include "Matrix.h" +using namespace std; + +void help_display(){ + cout<<"Help:"< input file (required) format .pdb"< normalized temperature (default =0.1)"< number of frames (default = -1)"< number of modes to consider (default = all)"< parameters file (DEFAULT = PARAMS_RIBOGM.DAT)"< PCA file (DEFAULT = none) format .pdb"< output file for fluctuations (DEFAULT = none)"< output file for correlations (DEFAULT = none)"< input file with shape data [data in third column] (DEFAULT = none)"< compute the Bhattacharya distance (DEFAULT = false)"< is a riboswitch? (do I need to take away the last res because it's an ADA?) (default=false)"< dump the ENM results on file (default=false)"< are there pre-constructed ENM files? (DEFAULT = false)"< (DEFAULT = intput file name)"< use C6 instead of C2 (default=false)"< bead_list; + bool out_has_name=false; + bool shape_has_name=false; + bool corr_has_name=false; + bool pca_has_name=false; + bool comp_bhatt=false; + bool ribo=false; + bool dump=false; + bool make_enm=true; + bool has_name=false; + bool use_C6=false; + // bead_list.push_back("C1'"); + // bead_list.push_back("C2"); + // bead_list.push_back("P"); + + cout<<" --- Program GenConf --- "<argc-1){cout<<"ERROR: empty parameter"<argc-1){cout<<"ERROR: empty parameter"<argc-1){cout<<"ERROR: empty parameter"<argc-1){cout<<"ERROR: empty parameter"<argc-1){cout<<"ERROR: empty parameter"<argc-1){cout<<"ERROR: empty parameter"<argc-1){cout<<"ERROR: empty parameter"<argc-1){cout<<"ERROR: empty parameter"<argc-1){cout<<"ERROR: empty parameter"<argc-1){cout<<"ERROR: empty parameter"< invalid parameter"< indexC2; + for (size_t i = 0; i < ref_struc.getSize(); ++i){ + if(ref_struc.getBead(i).getHet()) continue; + // if ((ref_struc.getBead(i).getAtomType()=="C2") || + // (ref_struc.getBead(i).getAtomType()=="C6") + if(use_C6){ + if (ref_struc.getBead(i).getAtomType()=="C6") + {is_C2[i]=1; indexC2.push_back(i);}// cout<<"eccolo: "< fluct; //vettore su cui salvero' le fluttuazioni C2-C2 + + if(NFRAMES<3){ //se il numero di frame e' 2 o meno non genero traiettorie e uso la formula "analitica" + if(NFRAMES==2) cout<<"Be careful! 2 frames doesn't make sense. I'm using the exact approximation"< nd(0.0,1); // definisco la distribuzione normale da cui estrarro i numeri gaussiani + boost::variate_generator > rnd(rng, nd); + + if((NTOP<1)||(NTOP>nmodes-6)) NTOP=nmodes-6; + vector sqrt_eval; + for (int m = 6; m < NTOP+6; ++m) { + double s=ENM.CovGetEigenval(nmodes-m-1); + if(s<0.0){cout<<"genConf *** WARNING: negative eigenvalue -> "< coords(indexC2.size()); + + int NFPRINT=NFRAMES/10; + if(NFPRINT<1) NFPRINT=1; + for (size_t f = 0; f < NFRAMES; ++f){ + if(f%NFPRINT==0) cout< width; + for (int m = 6; m < NTOP+6; ++m) { + //double s=ENM.CovGetEigenval(nmodes-m-1); + double s=sqrt_eval[m-6]; + width.push_back(rnd()*s); + } + // $ $ $ $ $ $ $ $ +#pragma omp parallel for + for (size_t i = 0; i < indexC2.size(); ++i){ + int i_bead=indexC2.at(i); + for (int m = 6; m < NTOP+6; ++m) { + coords[i] += ENM.CovGetEigenvec(nmodes-m-1)[i_bead]*width[m-6]; //sbagliavo qui... + } + } + // $ $ $ $ $ $ $ $ + + // # # # # # # # + for (size_t i = 0; i < indexC2.size(); ++i){ //### qui calcolo le distanze + double temp_dsq=coords[i].distanceSQ(coords[i+1]); + d_sq_mean[i]+=temp_dsq; + dist_mean[i]+=sqrt(temp_dsq); + } + + }//enddo frames + + for (size_t i = 0; i < n_C2; ++i){ + // for (size_t i = 0; i < indexC2.size()-2; ++i){ //NB: -2 per escludere ADA del riboswitch, altrimenti -1 !! + d_sq_mean[i]/=NFRAMES; + dist_mean[i]/=NFRAMES; + fluct.push_back(d_sq_mean[i]-(dist_mean[i]*dist_mean[i])); + } + }//endelse + + /// +++ OPEN the OUTPUT file +++ + if(out_has_name){ + ofstream fout; + fout.open(oname); + // ENM.setFast(false); + // ENM.computeCovar(); + + // ENM.Dump("prova"); + // ENM.dump_top_modes(10,5,"prova"); + // ### qui calcolo le distanze e le printo ### + // ENM.dumpDistFluc("prova"); + for (size_t i = 0; i < n_C2; ++i){ + // for (size_t i = 0; i < indexC2.size()-2; ++i){ //NB: -2 per escludere ADA del riboswitch, altrimenti -1 !! + fout< shape; + ifstream shapefile(shape_name); + char cdum[2]; + int idum; + double ddum; + double tmpshape; + /// while(shapefile >> cdum >> idum >> tmpshape >> cdum >> idum >>ddum >> cdum >> idum >>ddum >> cdum >> idum >>ddum >> cdum >> idum >>ddum){ + string line; + while(getline(shapefile,line)){ + std::istringstream iss(line); + iss >> cdum >> idum >> tmpshape; + if(!strncmp(cdum,"#",1)) continue; + shape.push_back(tmpshape); + } + cout<<"### corr ###"<0.0) + C+=1; + else + if((fluct[i]-fluct[j])*(shape[i]-shape[j])<0.0) + D+=1; + }//enddo j + }//enddo i + kend=(C-D)/C_n_2; + corr=(sum_XY-(sum_X*sum_Y)/ndata)/(sqrt((sum_X2-(sum_X*sum_X)/ndata)*(sum_Y2-(sum_Y*sum_Y)/ndata))); + cout< evec_ENM=ENM.CovGetEigenvec(m-6); + for (int n = 6; n < nmodes; ++n) { + double eval_PCA=PCA.getCovMat().GetEigenval(n); + vector evec_PCA=PCA.getCovMat().GetEigenvec(n); + double scal_prod=0; + for(int i=0; inu_frac*somma) break; + } + cout<<"NTOP="< +#include +#include +#include +#include "omp.h" +#include +#include + + +#include "Structure3d.h" +#include "PrincipalComp.h" +#include "ElasticNet.h" +//#include "Matrix.h" +using namespace std; + +#include "correlation.cc" + +void help_display(){ + cout<<"Help:"< input file (required) format .pdb"< number of modes to consider (default = all)"< parameters file (DEFAULT = PARAMS_RIBOGM.DAT)"< beads to consider (default = C1' C2 P )"< cutoff value in angstrom (default=10)"< multiple cutoff value in angstrom (default= single cutoff)"< cutoff value in angstrom (default=15), for protein residues"< beads to consider for output(default = same as beads)"< PCA file (DEFAULT = none) format .pdb"< output file for correlation and RWSIP (DEFAULT = none)"< compute the Bhattacharya distance (DEFAULT = false)"< is a riboswitch? (do I need to take away the last res because it's an ADA?) (default=false)"< dump the ENM results on file (default=false)"< are there pre-constructed ENM files? (DEFAULT = false)"< (DEFAULT = intput file name)"< compute force covariance"< compute covariance matrix"< bead_list; + vector outbead_list; + vector cutoff_list; + bool corr_has_name=false; + bool pca_has_name=false; + bool comp_bhatt=false; + bool ribo=false; + bool dump=false; + bool make_enm=true; + bool has_name=false; + bool compforce=false; + bool compcov=false; + bool par_has_name=false; + bool ENM_ERROR=false; + bool dump_nn=false; + bool multi_cutoff=false; + bool use_exp=false; + bead_list.push_back("C1'"); + bead_list.push_back("C2"); + bead_list.push_back("P"); + + cout<<" --- Program ENM_New --- "<argc-1){cout<<"ERROR: empty parameter"<argc-1){cout<<"ERROR: empty parameter"<argc-1){cout<<"ERROR: empty parameter"<argc-1){cout<<"ERROR: empty parameter"<argc-1){cout<<"ERROR: empty parameter"<argc-1){cout<<"ERROR: empty parameter"<argc-1){cout<<"ERROR: empty parameter"<argc-1){cout<<"ERROR: empty parameter"<argc-1){cout<<"ERROR: empty parameter"<argc-1){cout<<"ERROR: empty parameter"<argc-1){cout<<"ERROR: empty parameter"<argc-1){cout<<"ERROR: empty parameter"<argc-1){cout<<"ERROR: empty parameter"<argc-1){cout<<"ERROR: empty parameter"<argc-1){cout<<"ERROR: empty parameter"<argc-1){cout<<"ERROR: empty parameter"< invalid parameter"<0){ + cout<<"Out Beads to consider are: "; + for(int i=0; i res_list; + if(par_has_name){ //qui leggo l'eventuale file parametri (lista di OUTRES) + char stringa[1000],line[1000]; + ifstream file_in; + file_in.open(par_name); + if(!file_in.is_open()){ cout<<"Problem opening parameters file: "<>stringa)) continue; + if(!strncmp(stringa,"OUTRES:",6)){ + cout<> num){cout<0) {cout<<"setto outbead"<0) {cout<<"ERROR READING THE FILE"< input file (required) format .pdb"< PCA file (required) format .pdb"< parameters file (DEFAULT = PARAMS_RIBOGM.DAT)"< output file (DEFAULT = $input_file+_rwsip.dat)"< beads to consider (default = [ C1' C2 P ] )"<argc-1){cout<<"ERROR: empty parameter"<argc-1){cout<<"ERROR: empty parameter"<argc-1){cout<<"ERROR: empty parameter"<argc-1){cout<<"ERROR: empty parameter"< invalid parameter"<nu_frac*somma) break; + } + cout<<"NTOP="< + +#include "Matrix.h" +#include "Structure3d.h" +#include "ElasticNet.h" +#include "PrincipalComp.h" + +using namespace std; + + +void help_display(){ + cout<<"Help:"< input file (required) format .pdb"< PCA file (required) format .pdb"< parameters file (DEFAULT = PARAMS_RIBOGM.DAT)"< output file (DEFAULT = $input_file+_rwsip.dat)"< beads to consider (default = [ C1' C2 P ] )"<argc-1){cout<<"ERROR: empty parameter"<argc-1){cout<<"ERROR: empty parameter"<argc-1){cout<<"ERROR: empty parameter"<argc-1){cout<<"ERROR: empty parameter"< invalid parameter"< evec_ENM=ENM.getCovMatrix().GetEigenvec(m); + for (int n = 6; n < nmodes; ++n) { + double eval_PCA=PCA.getCovMat().GetEigenval(n); + vector evec_PCA=PCA.getCovMat().GetEigenvec(n); + double scal_prod=0; + for(int i=0; i +#include +#include +#include +#include +#include +#include + +//#include "my_malloc.h" //libraries to allocate vectors and matrices +//#include "io.h" //definitions of input and output subroutines +#include "Matrix.h" //definition of the class CMatrix +//#include "PrincipalComp.h" +#include "Structure3d.h" + +using namespace std; + +void help_display(){ + cout<<"Help:"< input file (required) format .pdb"< output file (required) format .pdb"< reference structure file format .pdb (default reference=first frame of input trajectory)"< beads to consider (default = ALL )"< fit only on the first/last 3 beads"< parameters file (DEFAULT = none)"< bead_list; + // list.push_back("C1'"); + // list.push_back("C2"); + // list.push_back("P"); + + + cout<<" --- Program FIT --- "<argc-1){cout<<"ERROR: empty parameter"<argc-1){cout<<"ERROR: empty parameter"<argc-1){cout<<"ERROR: empty parameter"<argc-1){cout<<"ERROR: empty parameter"<argc-1){cout<<"ERROR: empty parameter"<argc-1){cout<<"ERROR: empty parameter"<argc-1){cout<<"ERROR: empty parameter"<argc-1){cout<<"ERROR: empty parameter"<argc-1){cout<<"ERROR: empty parameter"< invalid parameter"< res_list; + if(par_has_name){ //qui leggo l'eventuale file parametri (lista di OUTRES) + char stringa[1000],line[1000]; + ifstream file_in; + file_in.open(par_name); + if(!file_in.is_open()){ cout<<"Problem opening parameters file: "<>stringa)) continue; + if(!strncmp(stringa,"OUTRES:",6)){ + cout<> num){cout< +#include +#include +#include + +#include "Structure3d.h" +#include "PrincipalComp.h" +//#include "Matrix.h" +using namespace std; + +void help_display(){ + cout<<"Help:"< input file (required) format .pdb"< normalized temperature (default =1)"< number of frames (default = 10000)"< number of modes to consider (default = all)"< beads to consider (default = [ C1' C2 P ] )"< bead_list; + bead_list.push_back("C1'"); + bead_list.push_back("C2"); + bead_list.push_back("P"); + + cout<<" --- Program GenConf --- "<argc-1){cout<<"ERROR: empty parameter"<argc-1){cout<<"ERROR: empty parameter"<argc-1){cout<<"ERROR: empty parameter"<argc-1){cout<<"ERROR: empty parameter"<argc-1){cout<<"ERROR: empty parameter"< invalid parameter"<0) nmodes=NTOP; + + cout<<"N. beads = "< coords(ref_struc.getSize()); + + // +++ OPEN the OUTPUT file +++ + char oname[1024]; + ofstream dumpfile; + sprintf(oname,"%s_trajectory.pdb",file_name); + dumpfile.open(oname); + // +++ +++ +++ +++ +++ +++ +++ + + int NFPRINT=NFRAMES/10; + if(NFPRINT<1) NFPRINT=1; + + for (size_t f = 0; f < NFRAMES; ++f){ + if(f%NFPRINT==0) cout<nmodes-6)) NTOP=nmodes-6; + for (int m = 0; m < NTOP-6; ++m) { + double s=nm.CovGetEigenval(nmodes-m-1); + // double s=nm.CovGetEigenval(m); + if(s<0.0001){ cout<<"genConf *** WARNING: negative eigenvalue -> "< nd(0.0,s); + + boost::variate_generator > rnd(rng, nd); + double width = rnd()/5.0*fattore; + for (size_t i = 0; i < ref_struc.getSize(); ++i){ + coords[i] += nm.CovGetEigenvec(nmodes-m-1)[i]*width; + } + } + + //printf("%ld\ncoordinates generated from normal modes\n",ref_struc.getSize()); + for (size_t i = 0; i < ref_struc.getSize(); ++i){ + char atomo[4]; + + new_struc.getBead(i).setCoordinates(coords[i]); + } + new_struc.dumpPDBFile(dumpfile); + + }//enddo frames + + return 0; +} diff --git a/io.h b/io.h new file mode 100644 index 0000000..8ff6e46 --- /dev/null +++ b/io.h @@ -0,0 +1,81 @@ +#ifndef IO_H_ +#define IO_H_ + +#include +#include + +/**********/ +static FILE *open_file_r (char *name){ + FILE *fp; + if ((fp = fopen (name, "r")) == NULL) + { + printf ("Could not open file %s.\n.Exiting.\n", name); + exit (1); + } + return (fp); +} +/**********/ +static FILE *open_file_w (char *name){ + FILE *fp; + if ((fp = fopen (name, "w")) == NULL) + { + printf ("Could not open file %s.\n.Exiting.\n", name); + exit (1); + } + return (fp); +} +/**********/ +static FILE *open_file_a (char *name){ + FILE *fp; + if ((fp = fopen (name, "a")) == NULL) + { + printf ("Could not open file %s.\n.Exiting.\n", name); + exit (1); + } + return (fp); +} +/**********/ +static void close_file (FILE * fp){ + fclose (fp); +} + +#endif /* IO_H_ */ + + +/**********/ + +//// ############ functions to read input from pdb ############ +//int get_mol_len_xyz(ifstream &ifile){ +// int mol_len=0; +// char line[200]; +// //find out the length of the molecule +// ifile>>mol_len; +// ifile.seekg(0, ios::beg); +// // ifile.getline(line,200); +// // cout<>idum; +// //cout<>chdum>>coordinates[3*i]>>coordinates[3*i+1]>>coordinates[3*i+2]; +// // cout< +#include +#include + +#include "kabsch.h" +#include "io.h" +#include "my_malloc.h" +#include "vectop.cc" +#include "myjacobi.cc" + +/* ============================================ */ +double rmsd_without_alignement(double str1[][3], double str2[][3],int length){ + /* ritorna l'rmsd tra str1 e str2. Non fa nessuna sovrapposizione ottimale */ + double rmsd, msd, v[3]; + int i,j; + + msd=0.0; + for(i=0 ; i < length; i++){ + for(j=0; j < 3; j++) { + v[j]= str1[i][j]-str2[i][j]; + msd+=v[j]*v[j]; + } + } + rmsd=sqrt(msd/length); + return(rmsd); +} + +/* ============================================ */ + +double drms(double str1[][3], double str2[][3],int length){ +/* ritorna l'rms delle distanze tra le N*(N-1) coppie di punti delle due strutture */ + + double drms, d1,d2; + int i,j; + + + drms=0.0; + for(i=0 ; i < length; i++){ + for(j=i+1 ; j < length; j++){ + d1 = dist_d(str1[i],str1[j],3); + d2 = dist_d(str2[i],str2[j],3); + drms += (d1-d2)*(d1-d2); + } + } + drms = drms*2/(length*(length-1)); + drms = sqrt(drms); + return(drms); +} + +/* ============================================ */ + +double best_rmsd(double str1[][3], double str2[][3],int length) +{ +/* ritorna il migliore rmsd (Kabsch) tra str1 e str2. Le coordinate di str1 e str2 vengono preservate */ + + double str3[2000][3], cm1[3],cm2[3], rmsd; + double u[3][3]; + + if (length > 2000) { + printf("aumenta la size vettore in best_rmsd\n"); + exit(1); + } + rmsd=optimal_alignment(str1,str2,str3,length,u,cm1,cm2); + return(rmsd); +} + +/* ============================================ */ +double kabsch(double str1[][3], double str2[][3], double str3[][3],int length) +{ +/* ritorna il migliore rmsd tra str1 e str2 senza allineare niente */ + + double rmsd; + double u[3][3], cm1[3],cm2[3]; + + rmsd=optimal_alignment(str1,str2,str3,length, u, cm1,cm2); + return(rmsd); + +} + +/* ============================================ */ + +double optimal_alignment(double str1[][3], double str2[][3], double str3[][3],int length, double u[][3], double *cm1, double *cm2) { + +/* ritorna il migliore rmsd tra str1 e str2 e mette in str3 la +prima struttura allineata sulla seconda, in u la matrice di rotazione */ + + void myjacobi (double a[][3], int n, double *d, double v[][3], int *nrot); + + int i, j, k, sign[3], order[3], nrot; + double e, e0; + double r[3][3], rt[3][3], temp, **x, **y; + double a[3][3], eval[3], evec[3][3]; + double eigenvalues[3], eigenvectors[3][3], b[3][3]; + + + x = d2t(length,3); + y = d2t(length,3); + + // zero_vec_d (cm1, 3); + // zero_vec_d (cm2, 3); + cm1=d1t(3); + cm2=d1t(3); + + + for (i = 0; i < length; i++) + { + for (j = 0; j < 3; j++) + { + cm1[j] += str1[i][j] / length; + cm2[j] += str2[i][j] / length; + + } + } +/* printf("cm1 %lf %lf %lf\n",cm1[0],cm1[1],cm1[2]); + printf("cm2 %lf %lf %lf\n",cm2[0],cm2[1],cm2[2]); +*/ + for (i = 0; i < length; i++) + { + for (j = 0; j < 3; j++) + { + x[i][j] = str1[i][j] - cm1[j]; + y[i][j] = str2[i][j] - cm2[j]; + } + } + + /* Mettiamo un pelo di noise per evitare allineamenti perfetti */ + + e0 = 0.0; + for (i = 0; i < length; i++) + { + e0 += 0.5 * norm_d (x[i], 3) * norm_d (x[i], 3); + e0 += 0.5 * norm_d (y[i], 3) * norm_d (y[i], 3); + } + + for (i = 0; i < 3; i++) + { + for (j = 0; j < 3; j++) + { + r[i][j] = 0.0; + for (k = 0; k < length; k++) + { + r[i][j] += y[k][i] * x[k][j]; + } + rt[j][i] = r[i][j]; + } + } + + + + for (i = 0; i < 3; i++) + { + for (j = 0; j < 3; j++) + { + a[i][j] = 0; + for (k = 0; k < 3; k++) + { + a[i][j] += rt[i][k] * r[k][j]; + } + } + } + + + myjacobi (a, 3, eval, evec, &nrot); +/* aggiungiamo delle piccole quantita' per rimuovere degenerazioni */ + eval[0]+=0.0000000000000001; + eval[1]+=0.00000000000000001; + eval[2]+=0.00000000000000002; + + if ((eval[0] < eval[1]) && (eval[0] < eval[2])) + { + order[0] = 1; + order[1] = 2; + order[2] = 0; + } + + if ((eval[1] < eval[0]) && (eval[1] < eval[2])) + { + order[0] = 0; + order[1] = 2; + order[2] = 1; + } + + if ((eval[2] < eval[0]) && (eval[2] < eval[1])) + { + order[0] = 0; + order[1] = 1; + order[2] = 2; + } + + + for (i = 0; i < 3; i++) + { + eigenvalues[i] = eval[order[i]]; + for (j = 0; j < 3; j++) + { + eigenvectors[i][j] = evec[j][order[i]]; + } + } + + + normalize_d (eigenvectors[0], 3); + normalize_d (eigenvectors[1], 3); + vecprod_d (eigenvectors[0], eigenvectors[1], eigenvectors[2]); + normalize_d (eigenvectors[2], 3); + + for (i = 0; i < 2; i++) + { + for (j = 0; j < 3; j++) + { + b[i][j] = 0; + for (k = 0; k < 3; k++) + { + b[i][j] += r[j][k] * eigenvectors[i][k]; + } + } + normalize_d (b[i], 3); + } + + + vecprod_d (b[0], b[1], b[2]); + normalize_d (b[2], 3); + + + temp = 0.0; + for (i = 0; i < 3; i++) + { + for (j = 0; j < 3; j++) + { + temp += b[2][i] * r[i][j] * eigenvectors[2][j]; + } + } + sign[2] = +1; + if (temp < 0) + sign[2] = -1; + sign[0]=sign[1]=1; + + e = e0 - sqrt (eigenvalues[0]) - sqrt (eigenvalues[1]) - sign[2] * sqrt (eigenvalues[2]); + + + +/********************/ + e =2.0 * e / length; + if (e <0.0) { + if (fabs(e) < 1.0e-3) { printf("Warning. In Kabsch alignment found slightly negative value of e (%e). Roundoff error? I will set it equal to zero.\n",e); e=0.0;}/* occasionally, when dealing with two practically identical configurations + the value of e may be slightly negative due to the small offsets and roundoff errors. + In this case we set it equal to zero. */ + else {fprintf(stderr,"ERROR. In Kabsch alignment found negative value of e: (%lf)\n",e); exit(1);} + } + e = sqrt (e); +/********************/ + + + for(i=0; i < 3; i++){ + for(j=0; j < 3; j++){ + u[i][j]=0.0; + for(k=0; k < 3; k++){ + u[i][j]+= b[k][i]*eigenvectors[k][j]; + + } +/* printf("%10.4lf ",u[i][j]); */ + + } +/* printf("\n"); */ + } + +/* allinea la 1 sulla 2 - compreso shift del centro di massa */ + for (i = 0; i < length; i++){ + for (j = 0; j < 3; j++){ + str3[i][j]=cm2[j]; + for (k = 0; k < 3; k++){ + str3[i][j]+=u[j][k]*x[i][k]; + } + } +/* printf("%8.3lf %8.3lf %8.3lf\n",str3[i][0],str3[i][1],str3[i][2]); + printf("AAAA %8.3lf %8.3lf %8.3lf\n",x[i][0],x[i][1],x[i][2]); +*/ + } + + free_d2t(x); + free_d2t(y); + return (e); + +} + + + +/* ============================================ */ + +double best_weighted_rmsd(double str1[][3], double str2[][3],double weight[],int length) +{ +/* ritorna il migliore rmsd tra str1 e str2 senza allineare niente */ + + double str3[2000][3], cm1[3],cm2[3], rmsd; + double u[3][3]; + + rmsd=optimal_weighted_alignment(str1,str2,str3,weight,length,u,cm1,cm2); + return(rmsd); +} + +/* ============================================ */ +double weighted_kabsch(double str1[][3], double str2[][3], double str3[][3],double weight[],int length) +{ +/* fa un allineamento ottimale pesato tra str1 e str2 e mette in str3 la prima struttura allineata sulla seconda. Ritorna l'rmsd pesato */ + + double rmsd; + double u[3][3], cm1[3],cm2[3]; + + rmsd=optimal_weighted_alignment(str1,str2,str3,weight,length, u, cm1,cm2); + return(rmsd); + +} + +/* ============================================ */ + + +double optimal_weighted_alignment(double str1[][3], double str2[][3], double str3[][3],double *weight,int length, double u[][3], double *cm1, double *cm2) { + +/*fa un allineamento ottimale pesato tra str1 e str2 e mette in str3 la +prima struttura allineata sulla seconda, in u la matrice di rotazione. Ritorna l'rmsd pesato */ + +/* Corretta il 11/6/2010. Non veniva considerato il peso per calcolare il centro di massa, cioe' l'origine attorno a cui ruotare le due strutture */ + + + void myjacobi (double a[][3], int n, double *d, double v[][3], int *nrot); + + int i, j, k, sign[3], order[3], nrot; + double e, e0; + double r[3][3], rt[3][3], temp, **x, **y; + double a[3][3], eval[3], evec[3][3]; + double eigenvalues[3], eigenvectors[3][3], b[3][3]; + double tot_weight; + + x = d2t(length,3); + y = d2t(length,3); + + // zero_vec_d (cm1, 3); + // zero_vec_d (cm2, 3); + cm1=d1t(3); + cm2=d1t(3); + + tot_weight=0.0; + for (i = 0; i < length; i++) tot_weight+=weight[i]; + + for (i = 0; i < length; i++) + { + for (j = 0; j < 3; j++) + { + cm1[j] += weight[i]*str1[i][j]/tot_weight; + cm2[j] += weight[i]*str2[i][j]/tot_weight; + + } + + } + + +/* printf("cm1 %lf %lf %lf\n",cm1[0],cm1[1],cm1[2]); + printf("cm2 %lf %lf %lf\n",cm2[0],cm2[1],cm2[2]); +*/ + for (i = 0; i < length; i++) + { + for (j = 0; j < 3; j++) + { + x[i][j] = str1[i][j] - cm1[j]; + y[i][j] = str2[i][j] - cm2[j]; + } + } + + + e0 = 0.0; + for (i = 0; i < length; i++) + { + e0 += 0.5 * weight[i]*norm_d (x[i], 3) * norm_d (x[i], 3); + e0 += 0.5 * weight[i]*norm_d (y[i], 3) * norm_d (y[i], 3); + } + + for (i = 0; i < 3; i++) + { + for (j = 0; j < 3; j++) + { + r[i][j] = 0.0; + for (k = 0; k < length; k++) + { + r[i][j] += weight[k]*y[k][i] * x[k][j]; + } + rt[j][i] = r[i][j]; + } + } + + + + for (i = 0; i < 3; i++) + { + for (j = 0; j < 3; j++) + { + a[i][j] = 0; + for (k = 0; k < 3; k++) + { + a[i][j] += rt[i][k] * r[k][j]; + } + } + } + + + myjacobi (a, 3, eval, evec, &nrot); +/* aggiungiamo delle piccole quantita' per rimuovere degenerazioni */ + eval[0]+=0.0000000000000001; + eval[1]+=0.00000000000000001; + eval[2]+=0.00000000000000002; + + if ((eval[0] < eval[1]) && (eval[0] < eval[2])) + { + order[0] = 1; + order[1] = 2; + order[2] = 0; + } + + if ((eval[1] < eval[0]) && (eval[1] < eval[2])) + { + order[0] = 0; + order[1] = 2; + order[2] = 1; + } + + if ((eval[2] < eval[0]) && (eval[2] < eval[1])) + { + order[0] = 0; + order[1] = 1; + order[2] = 2; + } + + + for (i = 0; i < 3; i++) + { + eigenvalues[i] = eval[order[i]]; + for (j = 0; j < 3; j++) + { + eigenvectors[i][j] = evec[j][order[i]]; + } + } + + + normalize_d (eigenvectors[0], 3); + normalize_d (eigenvectors[1], 3); + vecprod_d (eigenvectors[0], eigenvectors[1], eigenvectors[2]); + normalize_d (eigenvectors[2], 3); + + for (i = 0; i < 2; i++) + { + for (j = 0; j < 3; j++) + { + b[i][j] = 0; + for (k = 0; k < 3; k++) + { + b[i][j] += r[j][k] * eigenvectors[i][k]; + } + } + normalize_d (b[i], 3); + } + + + vecprod_d (b[0], b[1], b[2]); + normalize_d (b[2], 3); + + + temp = 0.0; + for (i = 0; i < 3; i++) + { + for (j = 0; j < 3; j++) + { + temp += b[2][i] * r[i][j] * eigenvectors[2][j]; + } + } + sign[2] = +1; + if (temp < 0) + sign[2] = -1; + sign[0]=sign[1]=1; + + e = e0 - sqrt (eigenvalues[0]) - sqrt (eigenvalues[1]) - sign[2] * sqrt (eigenvalues[2]); + + + +/********************/ + e =2.0 * e / tot_weight; + if (e <0.0) { + if (fabs(e) < 1.0e-3) { printf("Warning. In Kabsch alignment found slightly negative value of e (%e). Roundoff error? I will set it equal to zero.\n",e); e=0.0;}/* occasionally, when dealing with two practically identical configurations + the value of e may be slightly negative due to the small offsets and roundoff errors. + In this case we set it equal to zero. */ + else {fprintf(stderr,"ERROR. In Kabsch alignment found negative value of e: (%lf)\n",e); exit(1);} + } + e = sqrt (e); +/********************/ + + + + for(i=0; i < 3; i++){ + for(j=0; j < 3; j++){ + u[i][j]=0.0; + for(k=0; k < 3; k++){ + u[i][j]+= b[k][i]*eigenvectors[k][j]; + } +/* printf("%10.4lf ",u[i][j]); */ + + } +/* printf("\n"); */ + } + +/* allinea la 1 sulla 2 - compreso shift del centro di massa */ + for (i = 0; i < length; i++){ + for (j = 0; j < 3; j++){ + str3[i][j]=cm2[j]; + for (k = 0; k < 3; k++){ + str3[i][j]+=u[j][k]*x[i][k]; + } + } +/* printf("%8.3lf %8.3lf %8.3lf\n",str3[i][0],str3[i][1],str3[i][2]); + printf("AAAA %8.3lf %8.3lf %8.3lf\n",x[i][0],x[i][1],x[i][2]); +*/ + } + + free_d2t(x); + free_d2t(y); + return (e); + +} + + + diff --git a/kabsch.h b/kabsch.h new file mode 100644 index 0000000..a46644f --- /dev/null +++ b/kabsch.h @@ -0,0 +1,19 @@ +/* Written by C. Micheletti, michelet@sissa.it */ +/* Last revision June 2010 */ + + +double best_rmsd(double str1[][3], double str2[][3],int length); +/* migliore rmsd di allineamento tra 2 strutture */ +double kabsch(double str1[][3], double str2[][3], double str3[][3],int length); +/* da il migliore rmsd e mette in str3 l'allineamento della str1 sulla str2*/ +double optimal_alignment(double str1[][3], double str2[][3], double str3[][3],int length, double u[][3], double *cm1, double *cm2); +/* da il migliore rmsd e mette in str3 l'allineamento della str1 sulla str2,in u la matrice di rotazione e le traslazioni +del centro di massa da fare per portare la 1 sulla 2 */ + +double best_weighted_rmsd(double str1[][3], double str2[][3],double weight[],int length); +/* migliore rmsd pesato di allineamento tra 2 strutture */ +double weighted_kabsch(double str1[][3], double str2[][3], double str3[][3],double weight[],int length); +/* da il migliore rmsd pesato e mette in str3 l'allineamento della str1 sulla str2*/ +double optimal_weighted_alignment(double str1[][3], double str2[][3], double str3[][3],double weight[],int length, double u[][3], double *cm1, double *cm2); +/* da il migliore rmsd pesato e mette in str3 l'allineamento della str1 sulla str2,in u la matrice di rotazione e le traslazioni +del centro di massa da fare per portare la 1 sulla 2 */ diff --git a/lapack_matrix_routines_wrapper_v3.cc b/lapack_matrix_routines_wrapper_v3.cc new file mode 100644 index 0000000..641b90e --- /dev/null +++ b/lapack_matrix_routines_wrapper_v3.cc @@ -0,0 +1,252 @@ +extern "C" { +#include +} + + +//#include +#include +#include +//#include +//#include +#include + + +/* // copy this in you .h +extern "C" void dgetrf(int *,int *, double *, int *, int *, int *); +extern "C" void dgetri(int *,double *,int *,int *, double *, int*,int *); +extern "C" void dsyevd(char *,char *,int *, double *, int *,double *,double *, int *, int *,int *, int * ); +*/ + + +static int invert_symmetric_matrix_lapack( double **m, int dim_m, double **inv_m){ + +/* + M is a NON-SINGULAR symmetric matrix of size (dim_m * dim_m). + The program calculates the inverse and returns it in inv_m + + written by CM, 20/07/2008 + tested by CM, 21/07/2008 +*/ + + int c, i,j,k,l,lwork, info; + double *array_m; + double *work, tmp1, tmp2, worksize; + int *iwork; + int *isupps; + + // cout<<"cacca2"<= k 0 no */ + for( k = 0;k< dim_m; k++){ + /* Run the sum backwards to minimize roundoff errors in the sum of inverse eigenvalues */ + if (fabs(eigenvalues_m[k]) < tol) continue; + inv_m[i][j] += 1.0/eigenvalues_m[k]*eigenvectors[k][i]*eigenvectors[k][j]; + } + if (i!=j) inv_m[j][i]=inv_m[i][j]; /* symmetry */ + } + } + +} diff --git a/my_malloc.h b/my_malloc.h new file mode 100644 index 0000000..6f7744f --- /dev/null +++ b/my_malloc.h @@ -0,0 +1,870 @@ +#ifndef MY_MALLOC_H_ +#define MY_MALLOC_H_ + +#include +//#include +#include +#include +#include + +using namespace std; + +static void failed (string message) +{ + cout<<"\n FAILED *** "< *vd1t(int n1) +{ + vector *p; + if ((p = (vector *) malloc ((size_t) n1 * sizeof (vector))) == NULL) + failed ("vd1t: failed"); + return p; +} +//### + +static FILE **F2t (int n1) +{ + FILE **p; + if ((p = (FILE **) malloc ((size_t) n1 * sizeof (FILE *))) == NULL) + failed ("F2t: failed n1"); + return p; +} + +static FILE ***F3t (int n1, int n2) +{ + FILE ***p; + int i; + if ((p = (FILE ***) malloc ((size_t) n1 * sizeof (FILE **))) == NULL) + failed ("F3t: failed n1"); + if ((p[0] = (FILE **) malloc ((size_t) n1 * n2 * sizeof (FILE *))) == NULL) + failed ("F3t: failed n2"); + for (i = 0; i < n1 - 1; i++) + p[i + 1] = p[i] + n2; + return p; +} + +static char *c1t (int n1) +{ + char *p; + if ((p = (char *) malloc ((size_t) n1 * sizeof (char))) == NULL) + failed ("c1t: failed"); + return p; +} + +static char **c2t (int n1, int n2) +{ + char **p; + int i; + if ((p = (char **) malloc ((size_t) n1 * sizeof (char *))) == NULL) + failed ("c2t: failed n1"); + if ((p[0] = (char *) malloc ((size_t) n1 * n2 * sizeof (char))) == NULL) + failed ("c2t: failed n2"); + for (i = 0; i < n1 - 1; i++) + p[i + 1] = p[i] + n2; + return p; +} + +static char ***c3t (int n1, int n2, int n3) +{ + char ***p; + int i, j; + if ((p = (char ***) malloc ((size_t) n1 * sizeof (char **))) == NULL) + failed ("c3t: failed n1"); + if ((p[0] = (char **) malloc ((size_t) n1 * n2 * sizeof (char *))) == NULL) + failed ("c3t: failed n2"); + if ((p[0][0] = (char *) malloc ((size_t) n1 * n2 * n3 * sizeof (char))) == NULL) + failed ("c3t: failed n3"); + for (i = 0; i < n2 - 1; i++) + p[0][i + 1] = p[0][i] + n3; + for (i = 0; i < n1 - 1; i++) + { + p[i + 1] = p[i] + n2; + p[i + 1][0] = p[i][0] + n2 * n3; + for (j = 0; j < n2 - 1; j++) + p[i + 1][j + 1] = p[i + 1][j] + n3; + } + return p; +} + + +static char ****c4t (int n1, int n2, int n3, int n4) +{ + char ****p, *a; + int i, j, k; + if ((p = (char ****) malloc ((size_t) n1 * sizeof (char ***))) == NULL) + failed ("c4t: failed n1"); + if ((p[0] = (char ***) malloc ((size_t) n1 * n2 * sizeof (char **))) == NULL) + failed ("c4t: failed n2"); + if ((p[0][0] = (char **) malloc ((size_t) n1 * n2 * n3 * sizeof (char *))) == NULL) + failed ("c4t: failed n3"); + if ((p[0][0][0] = (char *) malloc ((size_t) n1 * n2 * n3 * n4 * sizeof (char ))) == NULL) + failed ("c4t: failed n4"); + + for (i = 0; i < n3 - 1; i++) + p[0][0][i + 1] = p[0][0][i] + n4; + + for (i = 0; i < n2 - 1; i++) + { + p[0][i + 1] = p[0][i] + n3; + p[0][i + 1][0] = p[0][i][0] + n3 * n4; + for (j = 0; j < n3 - 1; j++) + p[0][i + 1][j + 1] = p[0][i + 1][j] + n4; + } + + for (i = 0; i < n1 - 1; i++) + { + p[i + 1] = p[i] + n2; + p[i + 1][0] = p[i][0] + n2 * n3; + p[i + 1][0][0] = p[i][0][0] + n2 * n3 * n4; + + for (j = 0; j < n3 - 1; j++) + p[i + 1][0][j + 1] = p[i + 1][0][j] + n4; + + for (j = 0; j < n2 - 1; j++) + { + p[i + 1][j + 1] = p[i + 1][j] + n3; + p[i + 1][j + 1][0] = p[i + 1][j][0] + n3 * n4; + for (k = 0; k < n3 - 1; k++) + p[i + 1][j + 1][k + 1] = p[i + 1][j + 1][k] + n4; + } + } + + for (i = 0, a = p[0][0][0]; i < n1 * n2 * n3 * n4; i++) + *a++ = 0.0; + + return p; +} + + + + +static short *s1t (int n1) +{ + short *p, *a; + int i; + if ((p = (short *) malloc ((size_t) n1 * sizeof (short))) == NULL) + failed ("s1t: failed"); + for (i = 0, a = p; i < n1; i++) + *a++ = 0; + return p; +} + +static short **s2t (int n1, int n2) +{ + short **p, *a; + int i; + if ((p = (short **) malloc ((size_t) n1 * sizeof (short *))) == NULL) + failed ("s2t: failed n1"); + if ((p[0] = (short *) malloc ((size_t) n1 * n2 * sizeof (short))) == NULL) + failed ("s2t: failed n2"); + for (i = 0; i < n1 - 1; i++) + p[i + 1] = p[i] + n2; + for (i = 0, a = p[0]; i < n1 * n2; i++) + *a++ = 0; + return p; +} + +static short ***s3t (int n1, int n2, int n3) +{ + short ***p, *a; + int i, j; + if ((p = (short ***) malloc ((size_t) n1 * sizeof (short **))) == NULL) + failed ("s3t: failed n1"); + if ((p[0] = (short **) malloc ((size_t) n1 * n2 * sizeof (short *))) == NULL) + failed ("s3t: failed n2"); + if ((p[0][0] = (short *) malloc ((size_t) n1 * n2 * n3 * sizeof (short))) == NULL) + failed ("s3t: failed n3"); + for (i = 0; i < n2 - 1; i++) + p[0][i + 1] = p[0][i] + n3; + for (i = 0; i < n1 - 1; i++) + { + p[i + 1] = p[i] + n2; + p[i + 1][0] = p[i][0] + n2 * n3; + for (j = 0; j < n2 - 1; j++) + p[i + 1][j + 1] = p[i + 1][j] + n3; + } + for (i = 0, a = p[0][0]; i < n1 * n2 * n3; i++) + *a++ = 0; + return p; +} + +static short ****s4t (int n1, int n2, int n3, int n4) +{ + short ****p, *a; + int i, j, k; + if ((p = (short ****) malloc ((size_t) n1 * sizeof (short ***))) == NULL) + failed ("s4t: failed n1"); + if ((p[0] = (short ***) malloc ((size_t) n1 * n2 * sizeof (short **))) == NULL) + failed ("s4t: failed n2"); + if ((p[0][0] = (short **) malloc ((size_t) n1 * n2 * n3 * sizeof (short *))) == NULL) + failed ("s4t: failed n3"); + if ((p[0][0][0] = (short *) malloc ((size_t) n1 * n2 * n3 * n4 * sizeof (short ))) == NULL) + failed ("s4t: failed n4"); + + for (i = 0; i < n3 - 1; i++) + p[0][0][i + 1] = p[0][0][i] + n4; + + for (i = 0; i < n2 - 1; i++) + { + p[0][i + 1] = p[0][i] + n3; + p[0][i + 1][0] = p[0][i][0] + n3 * n4; + for (j = 0; j < n3 - 1; j++) + p[0][i + 1][j + 1] = p[0][i + 1][j] + n4; + } + + for (i = 0; i < n1 - 1; i++) + { + p[i + 1] = p[i] + n2; + p[i + 1][0] = p[i][0] + n2 * n3; + p[i + 1][0][0] = p[i][0][0] + n2 * n3 * n4; + + for (j = 0; j < n3 - 1; j++) + p[i + 1][0][j + 1] = p[i + 1][0][j] + n4; + + for (j = 0; j < n2 - 1; j++) + { + p[i + 1][j + 1] = p[i + 1][j] + n3; + p[i + 1][j + 1][0] = p[i + 1][j][0] + n3 * n4; + for (k = 0; k < n3 - 1; k++) + p[i + 1][j + 1][k + 1] = p[i + 1][j + 1][k] + n4; + } + } + + for (i = 0, a = p[0][0][0]; i < n1 * n2 * n3 * n4; i++) + *a++ = 0.0; + + return p; +} + + +static int *i1t (int n1){ + int *p, *a; + int i; + if ((p = (int *) malloc ((size_t) n1 * sizeof (int))) == NULL) + failed ("i1t: failed"); + for (i = 0, a = p; i < n1; i++) + *a++ = 0; + return p; +} + +static int **i2t (int n1, int n2){ + int **p, *a; + int i; + if ((p = (int **) malloc ((size_t) n1 * sizeof (int *))) == NULL) + failed ("i2t: failed n1"); + if ((p[0] = (int *) malloc ((size_t) n1 * n2 * sizeof (int))) == NULL) + failed ("i2t: failed n2"); + for (i = 0; i < n1 - 1; i++) + p[i + 1] = p[i] + n2; + for (i = 0, a = p[0]; i < n1 * n2; i++) + *a++ = 0; + return p; +} + +static int ***i3t (int n1, int n2, int n3){ + int ***p, *a; + int i, j; + if ((p = (int ***) malloc ((size_t) n1 * sizeof (int **))) == NULL) + failed ("i3t: failed n1"); + if ((p[0] = (int **) malloc ((size_t) n1 * n2 * sizeof (int *))) == NULL) + failed ("i3t: failed n2"); + if ((p[0][0] = (int *) malloc ((size_t) n1 * n2 * n3 * sizeof (int))) == NULL) + failed ("i3t: failed n3"); + for (i = 0; i < n2 - 1; i++) + p[0][i + 1] = p[0][i] + n3; + for (i = 0; i < n1 - 1; i++) + { + p[i + 1] = p[i] + n2; + p[i + 1][0] = p[i][0] + n2 * n3; + for (j = 0; j < n2 - 1; j++) + p[i + 1][j + 1] = p[i + 1][j] + n3; + } + for (i = 0, a = p[0][0]; i < n1 * n2 * n3; i++) + *a++ = 0; + return p; +} + +static int ****i4t (int n1, int n2, int n3, int n4){ + int ****p, *a; + int i, j, k; + if ((p = (int ****) malloc ((size_t) n1 * sizeof (int ***))) == NULL) + failed ("i4t: failed n1"); + if ((p[0] = (int ***) malloc ((size_t) n1 * n2 * sizeof (int **))) == NULL) + failed ("i4t: failed n2"); + if ((p[0][0] = (int **) malloc ((size_t) n1 * n2 * n3 * sizeof (int *))) == NULL) + failed ("i4t: failed n3"); + if ((p[0][0][0] = (int *) malloc ((size_t) n1 * n2 * n3 * n4 * sizeof (int ))) == NULL) + failed ("i4t: failed n4"); + + for (i = 0; i < n3 - 1; i++) + p[0][0][i + 1] = p[0][0][i] + n4; + + for (i = 0; i < n2 - 1; i++) + { + p[0][i + 1] = p[0][i] + n3; + p[0][i + 1][0] = p[0][i][0] + n3 * n4; + for (j = 0; j < n3 - 1; j++) + p[0][i + 1][j + 1] = p[0][i + 1][j] + n4; + } + + for (i = 0; i < n1 - 1; i++) + { + p[i + 1] = p[i] + n2; + p[i + 1][0] = p[i][0] + n2 * n3; + p[i + 1][0][0] = p[i][0][0] + n2 * n3 * n4; + + for (j = 0; j < n3 - 1; j++) + p[i + 1][0][j + 1] = p[i + 1][0][j] + n4; + + for (j = 0; j < n2 - 1; j++) + { + p[i + 1][j + 1] = p[i + 1][j] + n3; + p[i + 1][j + 1][0] = p[i + 1][j][0] + n3 * n4; + for (k = 0; k < n3 - 1; k++) + p[i + 1][j + 1][k + 1] = p[i + 1][j + 1][k] + n4; + } + } + + for (i = 0, a = p[0][0][0]; i < n1 * n2 * n3 * n4; i++) + *a++ = 0.0; + + return p; +} + + +static float *f1t (int n1){ + float *p, *a; + int i; + if ((p = (float *) malloc ((size_t) n1 * sizeof (float))) == NULL) + failed ("f1t: failed"); + for (i = 0, a = p; i < n1; i++) + *a++ = 0; + return p; +} + +static float **f2t (int n1, int n2){ + float **p, *a; + int i; + if ((p = (float **) malloc ((size_t) n1 * sizeof (float *))) == NULL) + failed ("f2t: failed n1"); + if ((p[0] = (float *) malloc ((size_t) n1 * n2 * sizeof (float))) == NULL) + failed ("f2t: failed n2"); + for (i = 0; i < n1 - 1; i++) + p[i + 1] = p[i] + n2; + for (i = 0, a = p[0]; i < n1 * n2; i++) + *a++ = 0; + return p; +} + +static float ***f3t (int n1, int n2, int n3){ + float ***p, *a; + int i, j; + if ((p = (float ***) malloc ((size_t) n1 * sizeof (float **))) == NULL) + failed ("f3t: failed n1"); + if ((p[0] = (float **) malloc ((size_t) n1 * n2 * sizeof (float *))) == NULL) + failed ("f3t: failed n2"); + if ((p[0][0] = (float *) malloc ((size_t) n1 * n2 * n3 * sizeof (float))) == NULL) + failed ("f3t: failed n3"); + for (i = 0; i < n2 - 1; i++) + p[0][i + 1] = p[0][i] + n3; + for (i = 0; i < n1 - 1; i++) + { + p[i + 1] = p[i] + n2; + p[i + 1][0] = p[i][0] + n2 * n3; + for (j = 0; j < n2 - 1; j++) + p[i + 1][j + 1] = p[i + 1][j] + n3; + } + for (i = 0, a = p[0][0]; i < n1 * n2 * n3; i++) + *a++ = 0; + return p; +} + +static float ****f4t (int n1, int n2, int n3, int n4){ + float ****p, *a; + int i, j, k; + if ((p = (float ****) malloc ((size_t) n1 * sizeof (float ***))) == NULL) + failed ("f4t: failed n1"); + if ((p[0] = (float ***) malloc ((size_t) n1 * n2 * sizeof (float **))) == NULL) + failed ("f4t: failed n2"); + if ((p[0][0] = (float **) malloc ((size_t) n1 * n2 * n3 * sizeof (float *))) == NULL) + failed ("f4t: failed n3"); + if ((p[0][0][0] = (float *) malloc ((size_t) n1 * n2 * n3 * n4 * sizeof (float))) == NULL) + failed ("f4t: failed n4"); + + for (i = 0; i < n3 - 1; i++) + p[0][0][i + 1] = p[0][0][i] + n4; + + for (i = 0; i < n2 - 1; i++) + { + p[0][i + 1] = p[0][i] + n3; + p[0][i + 1][0] = p[0][i][0] + n3 * n4; + for (j = 0; j < n3 - 1; j++) + p[0][i + 1][j + 1] = p[0][i + 1][j] + n4; + } + + for (i = 0; i < n1 - 1; i++) + { + p[i + 1] = p[i] + n2; + p[i + 1][0] = p[i][0] + n2 * n3; + p[i + 1][0][0] = p[i][0][0] + n2 * n3 * n4; + + for (j = 0; j < n3 - 1; j++) + p[i + 1][0][j + 1] = p[i + 1][0][j] + n4; + + for (j = 0; j < n2 - 1; j++) + { + p[i + 1][j + 1] = p[i + 1][j] + n3; + p[i + 1][j + 1][0] = p[i + 1][j][0] + n3 * n4; + for (k = 0; k < n3 - 1; k++) + p[i + 1][j + 1][k + 1] = p[i + 1][j + 1][k] + n4; + } + } + + for (i = 0, a = p[0][0][0]; i < n1 * n2 * n3 * n4; i++) + *a++ = 0.0; + + return p; +} + +static float *****f5t (int n1, int n2, int n3, int n4, int n5){ + float *****p, *a; + int i, j, k, l; + if ((p = (float *****) malloc ((size_t) n1 * sizeof (float ****))) == NULL) + failed ("f5t: failed n1"); + if ((p[0] = (float ****) malloc ((size_t) n1 * n2 * sizeof (float ***))) == NULL) + failed ("f5t: failed n2"); + if ((p[0][0] = + (float ***) malloc ((size_t) n1 * n2 * n3 * sizeof (float **))) == NULL) + failed ("f5t: failed n3"); + if ((p[0][0][0] = + (float **) malloc ((size_t) n1 * n2 * n3 * n4 * sizeof (float *))) == NULL) + failed ("f5t: failed n4"); + if ((p[0][0][0][0] = + (float *) malloc ((size_t) n1 * n2 * n3 * n4 * n5 * sizeof (float))) == NULL) + failed ("f5t: failed n5"); + + for (i = 0; i < n4 - 1; i++) + p[0][0][0][i + 1] = p[0][0][0][i] + n5; + + for (i = 0; i < n3 - 1; i++) + { + p[0][0][i + 1] = p[0][0][i] + n4; + p[0][0][i + 1][0] = p[0][0][i][0] + n4 * n5; + for (j = 0; j < n4 - 1; j++) + p[0][0][i + 1][j + 1] = p[0][0][i + 1][j] + n5; + } + + for (i = 0; i < n2 - 1; i++) + { + p[0][i + 1] = p[0][i] + n3; + p[0][i + 1][0] = p[0][i][0] + n3 * n4; + p[0][i + 1][0][0] = p[0][i][0][0] + n3 * n4 * n5; + + for (j = 0; j < n4 - 1; j++) + p[0][i + 1][0][j + 1] = p[0][i + 1][0][j] + n5; + + for (j = 0; j < n3 - 1; j++) + { + p[0][i + 1][j + 1] = p[0][i + 1][j] + n4; + p[0][i + 1][j + 1][0] = p[0][i + 1][j][0] + n4 * n5; + for (k = 0; k < n4 - 1; k++) + p[0][i + 1][j + 1][k + 1] = p[0][i + 1][j + 1][k] + n5; + } + } + + for (i = 0; i < n1 - 1; i++) + { + p[i + 1] = p[i] + n2; + p[i + 1][0] = p[i][0] + n2 * n3; + p[i + 1][0][0] = p[i][0][0] + n2 * n3 * n4; + p[i + 1][0][0][0] = p[i][0][0][0] + n2 * n3 * n4 * n5; + + for (j = 0; j < n4 - 1; j++) + p[i + 1][0][0][j + 1] = p[i + 1][0][0][j] + n5; + + for (j = 0; j < n3 - 1; j++) + { + p[i + 1][0][j + 1] = p[i + 1][0][j] + n4; + p[i + 1][0][j + 1][0] = p[i + 1][0][j][0] + n4 * n5; + for (k = 0; k < n4 - 1; k++) + p[i + 1][0][j + 1][k + 1] = p[i + 1][0][j + 1][k] + n5; + } + + for (j = 0; j < n2 - 1; j++) + { + p[i + 1][j + 1] = p[i + 1][j] + n3; + p[i + 1][j + 1][0] = p[i + 1][j][0] + n3 * n4; + p[i + 1][j + 1][0][0] = p[i + 1][j][0][0] + n3 * n4 * n5; + + for (k = 0; k < n3 - 1; k++) + p[i + 1][j + 1][0][k + 1] = p[i + 1][j + 1][0][k] + n5; + + for (k = 0; k < n3 - 1; k++) + { + p[i + 1][j + 1][k + 1] = p[i + 1][j + 1][k] + n4; + p[i + 1][j + 1][k + 1][0] = p[i + 1][j + 1][k][0] + n4 * n5; + for (l = 0; l < n4 - 1; l++) + p[i + 1][j + 1][k + 1][l + 1] = p[i + 1][j + 1][k + 1][l] + n5; + } + } + } + + for (i = 0, a = p[0][0][0][0]; i < n1 * n2 * n3 * n4 * n5; i++) + *a++ = 0.0; + + return p; +} + +static double *d1t (int n1){ + double *p, *a; + int i; + if ((p = (double *) malloc ((size_t) n1 * sizeof (double))) == NULL) + failed ("d1t: failed"); + for (i = 0, a = p; i < n1; i++) + *a++ = 0; + return p; +} + +static double **d2t (int n1, int n2){ + double **p, *a; + int i; + if ((p = (double **) malloc ((size_t) n1 * sizeof (double *))) == NULL) + failed ("d2t: failed n1"); + if ((p[0] = (double *) malloc ((size_t) n1 * n2 * sizeof (double))) == NULL) + failed ("d2t: failed n2"); + for (i = 0; i < n1 - 1; i++) + p[i + 1] = p[i] + n2; + for (i = 0, a = p[0]; i < n1 * n2; i++) + *a++ = 0; + return p; +} + +static double ***d3t (int n1, int n2, int n3){ + double ***p, *a; + int i, j; + if ((p = (double ***) malloc ((size_t) n1 * sizeof (double **))) == NULL) + failed ("d3t: failed n1"); + if ((p[0] = (double **) malloc ((size_t) n1 * n2 * sizeof (double *))) == NULL) + failed ("d3t: failed n2"); + if ((p[0][0] = (double *) malloc ((size_t) n1 * n2 * n3 * sizeof (double))) == NULL) + failed ("d3t: failed n3"); + for (i = 0; i < n2 - 1; i++) + p[0][i + 1] = p[0][i] + n3; + for (i = 0; i < n1 - 1; i++) + { + p[i + 1] = p[i] + n2; + p[i + 1][0] = p[i][0] + n2 * n3; + for (j = 0; j < n2 - 1; j++) + p[i + 1][j + 1] = p[i + 1][j] + n3; + } + for (i = 0, a = p[0][0]; i < n1 * n2 * n3; i++) + *a++ = 0; + return p; +} + + +static double ****d4t (int n1, int n2, int n3, int n4){ + double ****p, *a; + int i, j, k; + if ((p = (double ****) malloc ((size_t) n1 * sizeof (double ***))) == NULL) + failed ("d4t: failed n1"); + if ((p[0] = (double ***) malloc ((size_t) n1 * n2 * sizeof (double **))) == NULL) + failed ("d4t: failed n2"); + if ((p[0][0] = (double **) malloc ((size_t) n1 * n2 * n3 * sizeof (double *))) == NULL) + failed ("d4t: failed n3"); + if ((p[0][0][0] = (double *) malloc ((size_t) n1 * n2 * n3 * n4 * sizeof (double ))) == NULL) + failed ("d4t: failed n4"); + + for (i = 0; i < n3 - 1; i++) + p[0][0][i + 1] = p[0][0][i] + n4; + + for (i = 0; i < n2 - 1; i++){ + p[0][i + 1] = p[0][i] + n3; + p[0][i + 1][0] = p[0][i][0] + n3 * n4; + for (j = 0; j < n3 - 1; j++) + p[0][i + 1][j + 1] = p[0][i + 1][j] + n4; + } + + for (i = 0; i < n1 - 1; i++){ + p[i + 1] = p[i] + n2; + p[i + 1][0] = p[i][0] + n2 * n3; + p[i + 1][0][0] = p[i][0][0] + n2 * n3 * n4; + + for (j = 0; j < n3 - 1; j++) + p[i + 1][0][j + 1] = p[i + 1][0][j] + n4; + + for (j = 0; j < n2 - 1; j++){ + p[i + 1][j + 1] = p[i + 1][j] + n3; + p[i + 1][j + 1][0] = p[i + 1][j][0] + n3 * n4; + for (k = 0; k < n3 - 1; k++) + p[i + 1][j + 1][k + 1] = p[i + 1][j + 1][k] + n4; + } + } + + for (i = 0, a = p[0][0][0]; i < n1 * n2 * n3 * n4; i++) + *a++ = 0.0; + + return p; +} + +//###queste funzioni da qui in giu' a cosa mi servono??### +static void readeol (FILE * fp){ + char s; + while ((s = getc (fp)) != EOF) + if (s == '\n') + return; +} + + +static int myipow (int a, int n){ + int b = a; + while (--n > 0) + a *= b; + return a; +} + +static float myfpow (float a, int n){ + float b = a; + while (--n > 0) + a *= b; + return a; +} + +static double mydpow (double a, int n){ + double b = a; + while (--n > 0) + a *= b; + return a; +} + +static void pdarray (int n, int m, double **a){ + int i, j; + printf ("\n"); + for (i = 0; i < n; i++) + { + for (j = 0; j < m; j++) + printf ("%8.4f ", a[i][j]); + printf ("\n"); + } +} + +static void pdvector (int n, double *a){ + int i; + printf ("\n"); + for (i = 0; i < n; i++) + printf ("%8.4f ", a[i]); + printf ("\n"); +} + +static void pfarray (int n, int m, float **a){ + int i, j; + printf ("\n"); + for (i = 0; i < n; i++) + { + for (j = 0; j < m; j++) + printf ("%8.4f ", a[i][j]); + printf ("\n"); + } +} + +static void pfvector (int n, float *a){ + int i; + printf ("\n"); + for (i = 0; i < n; i++) + printf ("%8.4f ", a[i]); + printf ("\n"); +} + +static void piarray (int n, int m, int **a){ + int i, j; + printf ("\n"); + for (i = 0; i < n; i++) + { + for (j = 0; j < m; j++) + printf ("%8d ", a[i][j]); + printf ("\n"); + } +} + +static void pivector (int n, int *a){ + int i; + printf ("\n"); + for (i = 0; i < n; i++) + printf ("%8d ", a[i]); + printf ("\n"); +} + +static void zdarray (int n, int m, double **a){ + int i, j; + for (i = 0; i < n; i++) + for (j = 0; j < m; j++) + a[i][j] = 0.0; +} + +static void zdvector (int n, double *a){ + int i; + for (i = 0; i < n; i++) + a[i] = 0.0; +} + +static void zfarray (int n, int m, float **a){ + int i, j; + for (i = 0; i < n; i++) + for (j = 0; j < m; j++) + a[i][j] = 0.0; +} + +static void zfvector (int n, float *a){ + int i; + for (i = 0; i < n; i++) + a[i] = 0.0; +} + +static void ziarray (int n, int m, int **a){ + int i, j; + for (i = 0; i < n; i++) + for (j = 0; j < m; j++) + a[i][j] = 0; +} + +static void zivector (int n, int *a){ + int i; + for (i = 0; i < n; i++) + a[i] = 0; +} + + + +/*********************************/ + +static void free_c4t(char ****p){ + free(p[0][0][0]); + free(p[0][0]); + free(p[0]); + free(p); +} + +static void free_c3t(char ***p){ + free(p[0][0]); + free(p[0]); + free(p); +} + +static void free_c2t(char **p){ + free(p[0]); + free(p); +} + +static void free_c1t(char *p){ + free(p); +} + + +/*********************************/ + +static void free_s4t(short ****p){ + free(p[0][0][0]); + free(p[0][0]); + free(p[0]); + free(p); +} + +static void free_s3t(short ***p){ + free(p[0][0]); + free(p[0]); + free(p); +} + +static void free_s2t(short **p){ + free(p[0]); + free(p); +} + +static void free_s1t(short *p){ + free(p); +} +/*********************************/ + +static void free_i4t(int ****p){ + free(p[0][0][0]); + free(p[0][0]); + free(p[0]); + free(p); +} + +static void free_i3t(int ***p){ + free(p[0][0]); + free(p[0]); + free(p); +} + +static void free_i2t(int **p){ + free(p[0]); + free(p); +} + +static void free_i1t(int *p){ + free(p); +} + +/*********************************/ + +static void free_f5t(float *****p){ + free(p[0][0][0][0]); + free(p[0][0][0]); + free(p[0][0]); + free(p[0]); + free(p); +} + +static void free_f4t(float ****p){ + free(p[0][0][0]); + free(p[0][0]); + free(p[0]); + free(p); +} + +static void free_f3t(float ***p){ + free(p[0][0]); + free(p[0]); + free(p); +} + +static void free_f2t(float **p){ + free(p[0]); + free(p); +} + +static void free_f1t(float *p){ + free(p); +} + +/*********************************/ + +static void free_d4t(double ****p){ + free(p[0][0][0]); + free(p[0][0]); + free(p[0]); + free(p); +} + +static void free_d3t(double ***p){ + free(p[0][0]); + free(p[0]); + free(p); +} + +static void free_d2t(double **p){ + free(p[0]); + free(p); +} + +static void free_d1t(double *p){ + free(p); +} + +#endif /* MY_MALLOC_*/ diff --git a/myjacobi.cc b/myjacobi.cc new file mode 100644 index 0000000..f1dfdce --- /dev/null +++ b/myjacobi.cc @@ -0,0 +1,105 @@ +/* Written by C. Micheletti, michelet@sissa.it */ +/* Last revision January 2007 */ + +/* ============================================ */ + +void myjacobi (double a[][3], int n, double *d, double v[][3], int *nrot) +{ + +/* modificata 24/1/2007 per eliminare controlli di uguaglianza tra floats */ + + int j, iq, ip, i; + double tresh, theta, tau, t, sm, s, h, g, c; + double b[3], z[3]; + +#define ROTATE(a,i,j,k,l) g=a[i][j];h=a[k][l];a[i][j]=g-s*(h+g*tau);a[k][l]=h+s*(g-h*tau); + + for (ip = 0; ip <= (n - 1); ip++) + { + for (iq = 0; iq <= (n - 1); iq++) + v[ip][iq] = 0.0; + v[ip][ip] = 1.0; + } + for (ip = 0; ip <= (n - 1); ip++) + { + b[ip] = d[ip] = a[ip][ip]; + z[ip] = 0.0; + } + *nrot = 0; + for (i = 1; i <= 500; i++) + { + sm = 0.0; + for (ip = 0; ip <= n - 2; ip++) + { + for (iq = ip + 1; iq <= (n - 1); iq++) + sm += fabs (a[ip][iq]); + } + if (sm == 0.0) + { + return; + } + if (i < 4) + tresh = 0.2 * sm / (n * n); + else + tresh = 0.0; + for (ip = 0; ip <= n - 2; ip++) + { + for (iq = ip + 1; iq <= (n - 1); iq++) + { + g = 100.0 * fabs (a[ip][iq]); + if (i > 4 && (fabs ((fabs (d[ip]) + g) - fabs (d[ip])) < 1.0e-6) + && (fabs( (fabs(d[iq]) + g)- fabs (d[iq])) < 1.0e-6)) + a[ip][iq] = 0.0; + else if (fabs (a[ip][iq]) > tresh) + { + h = d[iq] - d[ip]; + if ( fabs((fabs (h) + g) - fabs (h)) < 1.0e-6) + t = (a[ip][iq]) / h; + else + { + theta = 0.5 * h / (a[ip][iq]); + t = 1.0 / (fabs (theta) + sqrt (1.0 + theta * theta)); + if (theta < 0.0) + t = -t; + } + c = 1.0 / sqrt (1 + t * t); + s = t * c; + tau = s / (1.0 + c); + h = t * a[ip][iq]; + z[ip] -= h; + z[iq] += h; + d[ip] -= h; + d[iq] += h; + a[ip][iq] = 0.0; + for (j = 0; j <= ip - 1; j++) + { + ROTATE (a, j, ip, j, iq) + } + for (j = ip + 1; j <= iq - 1; j++) + { + ROTATE (a, ip, j, j, iq) + } + for (j = iq + 1; j <= (n - 1); j++) + { + ROTATE (a, ip, j, iq, j) + } + for (j = 0; j <= (n - 1); j++) + { + ROTATE (v, j, ip, j, iq) + } + ++(*nrot); + } + } + } + for (ip = 0; ip <= (n - 1); ip++) + { + b[ip] += z[ip]; + d[ip] = b[ip]; + z[ip] = 0.0; + } + } + printf ("Too many iterations in routine JACOBI %lf",sm); + /* exit (1); */ +#undef ROTATE +} + diff --git a/pca.C b/pca.C new file mode 100644 index 0000000..48bc9c0 --- /dev/null +++ b/pca.C @@ -0,0 +1,251 @@ +/* Program written by Giovanni Pinamonti + PhD student at + Scuola Internazionale Superiori di Studi Avanzati, Trieste, Italy + begin june 11th 2013 */ + +#include +#include +#include +#include +#include +#include +#include + +//#include "my_malloc.h" //libraries to allocate vectors and matrices +//#include "io.h" //definitions of input and output subroutines +#include "Matrix.h" //definition of the class Matrix +#include "PrincipalComp.h" + + +using namespace std; + +void help_display(){ + cout<<"Help:"< input file (required)"< number of eigenvectors to write on file"< number of components on which project on"< beads to consider (default = ALL )"< skip the first/last residues (default = 0 )"< dump partial mean square fluctuations (default = false)"< compute effective interaction matrix? (default = false)"< dump the covariance matrix? (default = false)"< dump the fluctuations of distances? (default = false)"< maximum number of time steps to read (default = -1 i.e. all the trajectory)"< minimum number of time steps to start from (default = -1 i.e. all the trajectory)"< file to compute forces covariance (DEFAULT = none)"< (DEFAULT = intput file name)"< bead_list; + + + cout<<" --- Program PCA --- "<argc-1){cout<<"ERROR: empty parameter"<argc-1){cout<<"ERROR: empty parameter"<argc-1){cout<<"ERROR: empty parameter"<argc-1){cout<<"ERROR: empty parameter"<argc-1){cout<<"ERROR: empty parameter"<argc-1){cout<<"ERROR: empty parameter"<0) cout<<"stop after "<argc-1){cout<<"ERROR: empty parameter"<0) cout<<"stop after "<argc-1){cout<<"ERROR: empty parameter"<argc-1){cout<<"ERROR: empty parameter"< invalid parameter"<0){ + Structure3d temp_structure(file_name); + vector res_list; + int temp_size=temp_structure.getSizeRes(); + if( temp_size<2*nlast+1 ){ cout<<"ERROR: the structure in"<=nlast)&&(i0) cacca.set_ntmax(NTMAX); + if (NTMIN>0) cacca.set_ntmin(NTMIN); + cacca.readTrajectory(file_name); + // PrincipalComp cacca(file_name,lista); + //PrincipalComp cacca(file_name,bead_list); + // cacca.readTrajectory(file_name,bead_list); + // cacca(file_name); + cout<<"traj was read from -> "<0){ + cout< nor the names of its contributors may be used to + * endorse or promote products derived from this software without specific prior written + * permission. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS + * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT + * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A + * PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT + * HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, + * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT + * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, + * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY + * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE + * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + * + * Source: started anew. + * + * Change History: + * 2009/04/13 Started source + * 2010/03/28 Modified FastCalcRMSDAndRotation() to handle tiny qsqr + * If trying all rows of the adjoint still gives too small + * qsqr, then just return identity matrix. (DLT) + * 2010/06/30 Fixed prob in assigning A[9] = 0 in InnerProduct() + * invalid mem access + * 2011/02/21 Made CenterCoords use weights + * 2011/05/02 Finally changed CenterCoords declaration in qcprot.h + * Also changed some functions to static + * 2011/07/08 put in fabs() to fix taking sqrt of small neg numbers, fp error + * 2012/07/26 minor changes to comments and main.c, more info (v.1.4) + * + ******************************************************************************/ +#include "qcprot.h" + + +static double +InnerProduct(double *A, double **coords1, double **coords2, const int len, const double *weight) +{ + double x1, x2, y1, y2, z1, z2; + int i; + const double *fx1 = coords1[0], *fy1 = coords1[1], *fz1 = coords1[2]; + const double *fx2 = coords2[0], *fy2 = coords2[1], *fz2 = coords2[2]; + double G1 = 0.0, G2 = 0.0; + + A[0] = A[1] = A[2] = A[3] = A[4] = A[5] = A[6] = A[7] = A[8] = 0.0; + + if (weight != NULL) + { + for (i = 0; i < len; ++i) + { + x1 = weight[i] * fx1[i]; + y1 = weight[i] * fy1[i]; + z1 = weight[i] * fz1[i]; + + G1 += x1 * fx1[i] + y1 * fy1[i] + z1 * fz1[i]; + + x2 = fx2[i]; + y2 = fy2[i]; + z2 = fz2[i]; + + G2 += weight[i] * (x2 * x2 + y2 * y2 + z2 * z2); + + A[0] += (x1 * x2); + A[1] += (x1 * y2); + A[2] += (x1 * z2); + + A[3] += (y1 * x2); + A[4] += (y1 * y2); + A[5] += (y1 * z2); + + A[6] += (z1 * x2); + A[7] += (z1 * y2); + A[8] += (z1 * z2); + } + } + else + { + for (i = 0; i < len; ++i) + { + x1 = fx1[i]; + y1 = fy1[i]; + z1 = fz1[i]; + + G1 += x1 * x1 + y1 * y1 + z1 * z1; + + x2 = fx2[i]; + y2 = fy2[i]; + z2 = fz2[i]; + + G2 += (x2 * x2 + y2 * y2 + z2 * z2); + + A[0] += (x1 * x2); + A[1] += (x1 * y2); + A[2] += (x1 * z2); + + A[3] += (y1 * x2); + A[4] += (y1 * y2); + A[5] += (y1 * z2); + + A[6] += (z1 * x2); + A[7] += (z1 * y2); + A[8] += (z1 * z2); + } + } + + return (G1 + G2) * 0.5; +} + + +int +FastCalcRMSDAndRotation(Quaternion *q_max, double *A, double *rmsd, double E0, int len, double minScore) +{ + double Sxx, Sxy, Sxz, Syx, Syy, Syz, Szx, Szy, Szz; + double Szz2, Syy2, Sxx2, Sxy2, Syz2, Sxz2, Syx2, Szy2, Szx2, + SyzSzymSyySzz2, Sxx2Syy2Szz2Syz2Szy2, Sxy2Sxz2Syx2Szx2, + SxzpSzx, SyzpSzy, SxypSyx, SyzmSzy, + SxzmSzx, SxymSyx, SxxpSyy, SxxmSyy; + double C[4]; + int i; + double mxEigenV; + double oldg = 0.0; + double b, a, delta, rms, qsqr; + double q1, q2, q3, q4, normq; + double a11, a12, a13, a14, a21, a22, a23, a24; + double a31, a32, a33, a34, a41, a42, a43, a44; + double a2, x2, y2, z2; + double xy, az, zx, ay, yz, ax; + double a3344_4334, a3244_4234, a3243_4233, a3143_4133,a3144_4134, a3142_4132; + double evecprec = 1e-6; + double evalprec = 1e-11; + + Sxx = A[0]; Sxy = A[1]; Sxz = A[2]; + Syx = A[3]; Syy = A[4]; Syz = A[5]; + Szx = A[6]; Szy = A[7]; Szz = A[8]; + + Sxx2 = Sxx * Sxx; + Syy2 = Syy * Syy; + Szz2 = Szz * Szz; + + Sxy2 = Sxy * Sxy; + Syz2 = Syz * Syz; + Sxz2 = Sxz * Sxz; + + Syx2 = Syx * Syx; + Szy2 = Szy * Szy; + Szx2 = Szx * Szx; + + SyzSzymSyySzz2 = 2.0*(Syz*Szy - Syy*Szz); + Sxx2Syy2Szz2Syz2Szy2 = Syy2 + Szz2 - Sxx2 + Syz2 + Szy2; + + C[2] = -2.0 * (Sxx2 + Syy2 + Szz2 + Sxy2 + Syx2 + Sxz2 + Szx2 + Syz2 + Szy2); + C[1] = 8.0 * (Sxx*Syz*Szy + Syy*Szx*Sxz + Szz*Sxy*Syx - Sxx*Syy*Szz - Syz*Szx*Sxy - Szy*Syx*Sxz); + + SxzpSzx = Sxz + Szx; + SyzpSzy = Syz + Szy; + SxypSyx = Sxy + Syx; + SyzmSzy = Syz - Szy; + SxzmSzx = Sxz - Szx; + SxymSyx = Sxy - Syx; + SxxpSyy = Sxx + Syy; + SxxmSyy = Sxx - Syy; + Sxy2Sxz2Syx2Szx2 = Sxy2 + Sxz2 - Syx2 - Szx2; + + C[0] = Sxy2Sxz2Syx2Szx2 * Sxy2Sxz2Syx2Szx2 + + (Sxx2Syy2Szz2Syz2Szy2 + SyzSzymSyySzz2) * (Sxx2Syy2Szz2Syz2Szy2 - SyzSzymSyySzz2) + + (-(SxzpSzx)*(SyzmSzy)+(SxymSyx)*(SxxmSyy-Szz)) * (-(SxzmSzx)*(SyzpSzy)+(SxymSyx)*(SxxmSyy+Szz)) + + (-(SxzpSzx)*(SyzpSzy)-(SxypSyx)*(SxxpSyy-Szz)) * (-(SxzmSzx)*(SyzmSzy)-(SxypSyx)*(SxxpSyy+Szz)) + + (+(SxypSyx)*(SyzpSzy)+(SxzpSzx)*(SxxmSyy+Szz)) * (-(SxymSyx)*(SyzmSzy)+(SxzpSzx)*(SxxpSyy+Szz)) + + (+(SxypSyx)*(SyzmSzy)+(SxzmSzx)*(SxxmSyy-Szz)) * (-(SxymSyx)*(SyzpSzy)+(SxzmSzx)*(SxxpSyy-Szz)); + + /* Newton-Raphson */ + mxEigenV = E0; + for (i = 0; i < 50; ++i) + { + oldg = mxEigenV; + x2 = mxEigenV*mxEigenV; + b = (x2 + C[2])*mxEigenV; + a = b + C[1]; + delta = ((a*mxEigenV + C[0])/(2.0*x2*mxEigenV + b + a)); + mxEigenV -= delta; + /* printf("\n diff[%3d]: %16g %16g %16g", i, mxEigenV - oldg, evalprec*mxEigenV, mxEigenV); */ + if (fabs(mxEigenV - oldg) < fabs(evalprec*mxEigenV)) + break; + } + + if (i == 50) + fprintf(stderr,"\nMore than %d iterations needed!\n", i); + + /* the fabs() is to guard against extremely small, but *negative* numbers due to floating point error */ + rms = sqrt(fabs(2.0 * (E0 - mxEigenV)/len)); + (*rmsd) = rms; + /* printf("\n\n %16g %16g %16g \n", rms, E0, 2.0 * (E0 - mxEigenV)/len); */ + + if (minScore > 0) + if (rms < minScore) + return (-1); // Don't bother with rotation. + + a11 = SxxpSyy + Szz-mxEigenV; a12 = SyzmSzy; a13 = - SxzmSzx; a14 = SxymSyx; + a21 = SyzmSzy; a22 = SxxmSyy - Szz-mxEigenV; a23 = SxypSyx; a24= SxzpSzx; + a31 = a13; a32 = a23; a33 = Syy-Sxx-Szz - mxEigenV; a34 = SyzpSzy; + a41 = a14; a42 = a24; a43 = a34; a44 = Szz - SxxpSyy - mxEigenV; + a3344_4334 = a33 * a44 - a43 * a34; a3244_4234 = a32 * a44-a42*a34; + a3243_4233 = a32 * a43 - a42 * a33; a3143_4133 = a31 * a43-a41*a33; + a3144_4134 = a31 * a44 - a41 * a34; a3142_4132 = a31 * a42-a41*a32; + + q1 = a22*a3344_4334-a23*a3244_4234+a24*a3243_4233;// ### qui potrei + q2 = -a21*a3344_4334+a23*a3144_4134-a24*a3143_4133;// ### fare tutto + q3 = a21*a3244_4234-a22*a3144_4134+a24*a3142_4132;// ### con un + q4 = -a21*a3243_4233+a22*a3143_4133-a23*a3142_4132;// ### quaternione + + qsqr = q1 * q1 + q2 * q2 + q3 * q3 + q4 * q4; // ### pure qui + +/* The following code tries to calculate another column in the adjoint matrix when the norm of the + current column is too small. + Usually this block will never be activated. To be absolutely safe this should be + uncommented, but it is most likely unnecessary. +*/ + if (qsqr < evecprec) + { + q1 = a12*a3344_4334 - a13*a3244_4234 + a14*a3243_4233; + q2 = -a11*a3344_4334 + a13*a3144_4134 - a14*a3143_4133; + q3 = a11*a3244_4234 - a12*a3144_4134 + a14*a3142_4132; + q4 = -a11*a3243_4233 + a12*a3143_4133 - a13*a3142_4132; + qsqr = q1*q1 + q2 *q2 + q3*q3+q4*q4; + + if (qsqr < evecprec) + { + double a1324_1423 = a13 * a24 - a14 * a23, a1224_1422 = a12 * a24 - a14 * a22; + double a1223_1322 = a12 * a23 - a13 * a22, a1124_1421 = a11 * a24 - a14 * a21; + double a1123_1321 = a11 * a23 - a13 * a21, a1122_1221 = a11 * a22 - a12 * a21; + + q1 = a42 * a1324_1423 - a43 * a1224_1422 + a44 * a1223_1322; + q2 = -a41 * a1324_1423 + a43 * a1124_1421 - a44 * a1123_1321; + q3 = a41 * a1224_1422 - a42 * a1124_1421 + a44 * a1122_1221; + q4 = -a41 * a1223_1322 + a42 * a1123_1321 - a43 * a1122_1221; + qsqr = q1*q1 + q2 *q2 + q3*q3+q4*q4; + + if (qsqr < evecprec) + { + q1 = a32 * a1324_1423 - a33 * a1224_1422 + a34 * a1223_1322; + q2 = -a31 * a1324_1423 + a33 * a1124_1421 - a34 * a1123_1321; + q3 = a31 * a1224_1422 - a32 * a1124_1421 + a34 * a1122_1221; + q4 = -a31 * a1223_1322 + a32 * a1123_1321 - a33 * a1122_1221; + qsqr = q1*q1 + q2 *q2 + q3*q3 + q4*q4; + + if (qsqr < evecprec) + { + /* if qsqr is still too small, return the identity matrix. */ + // rot[0] = rot[4] = rot[8] = 1.0; + // rot[1] = rot[2] = rot[3] = rot[5] = rot[6] = rot[7] = 0.0; + + // return(0); + cerr<(q4,q1,q2,q3); + + /* QUI C'ERA UN MALEDETTO SEGNO MENO CHE NON TORNAVA!!! */ + (*q_max)=Quaternion(-q1,q2,q3,q4); + // (*q_max)=Quaternion(q1,q2,q3,q4); + + // normq = sqrt(qsqr); + // q1 /= normq; + // q2 /= normq; + // q3 /= normq; + // q4 /= normq; + + // a2 = q1 * q1; + // x2 = q2 * q2; + // y2 = q3 * q3; + // z2 = q4 * q4; + + // xy = q2 * q3; + // az = q1 * q4; + // zx = q4 * q2; + // ay = q1 * q3; + // yz = q3 * q4; + // ax = q1 * q2; + + // rot[0] = a2 + x2 - y2 - z2; + // rot[1] = 2 * (xy + az); + // rot[2] = 2 * (zx - ay); + // rot[3] = 2 * (xy - az); + // rot[4] = a2 - x2 + y2 - z2; + // rot[5] = 2 * (yz + ax); + // rot[6] = 2 * (zx + ay); + // rot[7] = 2 * (yz - ax); + // rot[8] = a2 - x2 - y2 + z2; + + return (1); +} + + +static Vector3d CenterCoords(double **coords, const int len, const double *weight){ + Vector3d com(0.,0.,0.); + double wsum; + double *x = coords[0], *y = coords[1], *z = coords[2]; + + if (weight != NULL){ + wsum = 0.0; + for (int i = 0; i < len; ++i){ + com.X += weight[i]*x[i]; + com.Y += weight[i]*y[i]; + com.Z += weight[i]*z[i]; + wsum += weight[i]; + } + com /= wsum; + } + else{ + for (int i = 0; i < len; ++i){ + com.X+= x[i]; + com.Y+= y[i]; + com.Z+= z[i]; + } + com /= len; + } + + for (int i = 0; i < len; ++i){ + x[i] -= com.X; + y[i] -= com.Y; + z[i] -= com.Z; + } + // cout< *q_max, const double *weight){ + double A[9], rmsd; + // /* center the structures -- if precentered you can omit this step */ + // Vector3d com1= CenterCoords(coords1, len, weight); + // Vector3d com2= CenterCoords(coords2, len, weight); + + /* calculate the (weighted) inner product of two structures */ + double E0 = InnerProduct(A, coords1, coords2, len, weight); + + /* calculate the RMSD & rotational matrix */ + FastCalcRMSDAndRotation(q_max, A, &rmsd, E0, len, -1); + + return rmsd; +} + + + + + +/* Function created by Giovanni Pinamonti + * it uses the QCP method to align coords2 (3Xlen) on coords1 (3Xlen) + * using the given weights + * and returns the best fit-RMSD; + */ +double QCPAlignment(double **coords1, double **coords2,const int len, const double *weight){ + + /* center the structures */ + Vector3d com1= CenterCoords(coords1,len,weight); + Vector3d com2= CenterCoords(coords2,len,weight); + + /* compute the RMSD and the rotational quaternion */ + Quaternion q_max; + double rmsd = CalcRMSDRotationalMatrix(coords1,coords2,len,&q_max,weight); + // l'rmsd lo fa giusto... + + /* apply the rotation to coords2 */ + for(int i = 0; i nor the names of its contributors may be used to + * endorse or promote products derived from this software without specific prior written + * permission. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS + * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT + * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A + * PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT + * HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, + * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT + * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, + * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY + * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE + * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + * + * Source: started anew. + * + * Change History: + * 2009/04/13 Started source + * + ******************************************************************************/ + +#include +#include +#include + +#include "quaternion.h" +#include "Vector3d.h" + +/* Calculate the RMSD & rotational matrix. + + Input: + coords1 -- reference structure + coords2 -- candidate structure + len -- the size of the system + weight -- the weight array of size len; set to NULL if not needed + Output: + rot[9] -- rotation matrix + Return: + RMSD value + +*/ +double CalcRMSDRotationalMatrix(double **coords1, double **coords2, const int len, Quaternion *q_max, const double *weight); + +/* Calculate the inner product of two structures. + If weight array is not NULL, calculate the weighted inner product. + + Input: + coords1 -- reference structure + coords2 -- candidate structure + len -- the size of the system + weight -- the weight array of size len: set to NULL if not needed + Output: + A[9] -- the inner product matrix + Return: + (G1 + G2) * 0.5; used as E0 in function 'FastCalcRMSDAndRotation' + + Warning: + 1. You MUST center the structures, coords1 and coords2, before calling this function. + + 2. Please note how the structure coordinates are stored in the double **coords + arrays. They are 3xN arrays, not Nx3 arrays as is also commonly + used (where the x, y, z axes are interleaved). The difference is + something like this for storage of a structure with 8 atoms: + + Nx3: xyzxyzxyzxyzxyzxyzxyzxyz + 3xN: xxxxxxxxyyyyyyyyzzzzzzzz + + The functions can be easily modified, however, to accomodate any + data format preference. I chose this format because it is readily + used in vectorized functions (SIMD, Altivec, MMX, SSE2, etc.). +*/ +//double InnerProduct(double *A, double **coords1, double **coords2, const int len, const double *weight); + +/* Calculate the RMSD, and/or the optimal rotation matrix. + + Input: + A[9] -- the inner product of two structures + E0 -- (G1 + G2) * 0.5 + len -- the size of the system + minScore-- if( minScore > 0 && rmsd < minScore) then calculate only the rmsd; + otherwise, calculate both the RMSD & the rotation matrix + Output: + rot[9] -- the rotation matrix in the order of xx, xy, xz, yx, yy, yz, zx, zy, zz + rmsd -- the RMSD value + Return: + only the rmsd was calculated if < 0 + both the RMSD & rotational matrix calculated if > 0 +*/ +int FastCalcRMSDAndRotation(Quaternion*q_max, double *A, double *rmsd, double E0, int len, double minScore); + +/* Center the coordinates. + + Warning: + If you are doing a full superposition (the usual least squares way), + you MUST center each structure first. That is, you must translate + each structure so that its centroid is at the origin. + You can use CenterCoords() for this. +*/ +//void CenterCoords(double **coords, const int len); + + + + +/* ### Questa l'ho fatta io ### */ +double QCPAlignment(double **,double**,const int,const double *); diff --git a/quaternion.c++ b/quaternion.c++ new file mode 100644 index 0000000..0654f2e --- /dev/null +++ b/quaternion.c++ @@ -0,0 +1,342 @@ +//**************************************************** +//* quaternion.c++ * +//* * +//* Implementaion for a generalized quaternion class * +//* * +//* Written 1.25.00 by Angela Bennett * +//**************************************************** + + +#include "quaternion.h" + +//Quaternion +// -default constructor +// -creates a new quaternion with all parts equal to zero +template +Quaternion<_Tp>::Quaternion(void) +{ + x = 0; + y = 0; + z = 0; + w = 0; +} + + +//Quaternion +// -constructor +// -parametes : x, y, z, w elements of the quaternion +// -creates a new quaternion based on the elements passed in +template +Quaternion<_Tp>::Quaternion(_Tp wi, _Tp xi, _Tp yi, _Tp zi) +{ + w = wi; + x = xi; + y = yi; + z = zi; +} + + +//Quaternion +// -constructor +// -parameters : vector/array of four elements +// -creates a new quaternion based on the elements passed in +template +Quaternion<_Tp>::Quaternion(_Tp v[4]) +{ + w = v[0]; + x = v[1]; + y = v[2]; + z = v[3]; +} + + +//Quaternion +// -copy constructor +// -parameters : const quaternion q +// -creates a new quaternion based on the quaternion passed in +template +Quaternion<_Tp>::Quaternion(const Quaternion<_Tp>& q) +{ + w = q.w; + x = q.x; + y = q.y; + z = q.z; +} + +#ifdef SHOEMAKE +//Quaternion +// -constructor +// -parameters : yaw, pitch, and roll of an Euler angle +// -creates a new quaternion based on the Euler elements passed in +// -used with Shoemakes code +template +Quaternion<_Tp>::Quaternion(_Tp e[3], int order) +{ + EulerAngles ea; + ea.x = e[0]; + ea.y = e[1]; + ea.z = e[2]; + ea.w = order; + + Quat q = Eul_ToQuat(ea); + + x = q.x; + y = q.y; + z = q.z; + w = q.w; +} +#endif + +//~Quaternion +// -destructor +// -deleted dynamically allocated memory +template +Quaternion<_Tp>::~Quaternion() +{ +} + + +//operator= +// -parameters : q1 - Quaternion object +// -return value : Quaternion +// -when called on quaternion q2 sets q2 to be an object of q3 +template +Quaternion<_Tp> Quaternion<_Tp>::operator = (const Quaternion<_Tp>& q) +{ + w = q.w; + x = q.x; + y = q.y; + z = q.z; + + return (*this); +} + +//operator+ +// -parameters : q1 - Quaternion object +// -return value : Quaternion +// -when called on quaternion q2 adds q1 + q2 and returns the sum in a new quaternion +template +Quaternion<_Tp> Quaternion<_Tp>::operator + (const Quaternion<_Tp>& q) +{ + return Quaternion(w+q.w, x+q.x, y+q.y, z+q.z); +} + +//operator- +// -parameters : q1- Quaternion object +// -return values : Quaternion +// -when called on q1 subtracts q1 - q2 and returns the difference as a new quaternion +template +Quaternion<_Tp> Quaternion<_Tp>::operator - (const Quaternion<_Tp>& q) +{ + return Quaternion(w-q.w, x-q.x, y-q.y, z-q.z); +} + + +//operator* +// -parameters : q1 - Quaternion object +// -return values : Quaternion +// -when called on a quaternion q2, multiplies q2 *q1 and returns the product in a new quaternion +template +Quaternion<_Tp> Quaternion<_Tp>::operator * (const Quaternion<_Tp>& q) +{ + return Quaternion( + w*q.w - x*q.x - y*q.y - z*q.z, + w*q.x + x*q.w + y*q.z - z*q.y, + w*q.y + y*q.w + z*q.x - x*q.z, + w*q.z + z*q.w + x*q.y - y*q.x); +} + +//operator/ +// -parameters : q1 and q2- Quaternion objects +// -return values : Quaternion +// -divide q1 by q2 and returns the quotient q1 +template +Quaternion<_Tp> Quaternion<_Tp>::operator / (Quaternion<_Tp>& q) +{ + return ((*this) * (q.inverse())); +} + + +//operator+= +// -parameters : q1- Quaternion object +// -return values : Quaternion +// -when called on quaternion q3, adds q1 and q3 and returns the sum as q3 +template +Quaternion<_Tp>& Quaternion<_Tp>::operator += (const Quaternion<_Tp>& q) +{ + w += q.w; + x += q.x; + y += q.y; + z += q.z; + + return (*this); +} + + +//operator-= +// -parameters : q1- Quaternion object +// -return values : Quaternion +// -when called on quaternion q3, subtracts q1 from q3 and returns the difference as q3 +template +Quaternion<_Tp>& Quaternion<_Tp>::operator -= (const Quaternion<_Tp>& q) +{ + w -= q.w; + x -= q.x; + y -= q.y; + z -= q.z; + + return (*this); +} + + +//operator*= +// -parameters : q1- Quaternion object +// -return values : Quaternion +// -when called on quaternion q3, multiplies q3 by q1 and returns the product as q3 +template +Quaternion<_Tp>& Quaternion<_Tp>::operator *= (const Quaternion<_Tp>& q) +{ + _Tp w_val = w*q.w - x*q.x - y*q.y - z*q.z; + _Tp x_val = w*q.x + x*q.w + y*q.z - z*q.y; + _Tp y_val = w*q.y + y*q.w + z*q.x - x*q.z; + _Tp z_val = w*q.z + z*q.w + x*q.y - y*q.x; + + w = w_val; + x = x_val; + y = y_val; + z = z_val; + + return (*this); +} + + +//operator/= +// -parameters : q1- Quaternion object +// -return values : quaternion +// -when called on quaternion q3, divides q3 by q1 and returns the quotient as q3 +template +Quaternion<_Tp>& Quaternion<_Tp>::operator /= (Quaternion<_Tp>& q) +{ + (*this) = (*this)*q.inverse(); + return (*this); +} + + +//operator!= +// -parameters : q1 and q2- Quaternion objects +// -return value : bool +// -determines if q1 and q2 are not equal +template +bool Quaternion<_Tp>::operator != (const Quaternion<_Tp>& q) +{ + return (w!=q.w || x!=q.x || y!=q.y || z!=q.z) ? true : false; +} + +//operator== +// -parameters : q1 and q2- Quaternion objects +// -return value : bool +// -determines if q1 and q2 are equal +template +bool Quaternion<_Tp>::operator == (const Quaternion<_Tp>& q) +{ + return (w==q.w && x==q.x && y==q.y && z==q.z) ? true : false; +} + +//norm +// -parameters : none +// -return value : _Tp +// -when called on a quaternion object q, returns the norm of q +template +_Tp Quaternion<_Tp>::norm() +{ + return (w*w + x*x + y*y + z*z); +} + +//magnitude +// -parameters : none +// -return value : _Tp +// -when called on a quaternion object q, returns the magnitude q +template +_Tp Quaternion<_Tp>::magnitude() +{ + return sqrt(norm()); +} + +//scale +// -parameters : s- a value to scale q1 by +// -return value: quaternion +// -returns the original quaternion with each part, w,x,y,z, multiplied by some scalar s +template +Quaternion<_Tp> Quaternion<_Tp>::scale(_Tp s) +{ + return Quaternion(w*s, x*s, y*s, z*s); +} + +// -parameters : none +// -return value : quaternion +// -when called on a quaternion object q, returns the inverse of q +template +Quaternion<_Tp> Quaternion<_Tp>::inverse() +{ + return conjugate().scale(1/norm()); +} + +//conjugate +// -parameters : none +// -return value : quaternion +// -when called on a quaternion object q, returns the conjugate of q +template +Quaternion<_Tp> Quaternion<_Tp>::conjugate() +{ + return Quaternion(w, -x, -y, -z); +} + +//UnitQuaternion +// -parameters : none +// -return value : quaternion +// -when called on quaterion q, takes q and returns the unit quaternion of q +template +Quaternion<_Tp> Quaternion<_Tp>::UnitQuaternion() +{ + return (*this).scale(1/(*this).magnitude()); +} + +// -parameters : vector of type _Tp +// -return value : void +// -when given a 3D vector, v, rotates v by this quaternion +template +void Quaternion<_Tp>::QuatRotation(_Tp v[3]) +{ + Quaternion <_Tp> qv(0, v[0], v[1], v[2]); + Quaternion <_Tp> qm = (*this) * qv * (*this).inverse(); + /* QUESTION: (*this) has to be normalized?? */ + v[0] = qm.x; + v[1] = qm.y; + v[2] = qm.z; +} + +#ifdef SHOEMAKE +// -parameters : integer order- which will specify the order of the rotation, q- quaternion +// -return value : Euler angle +// - +template +void Quaternion<_Tp>::toEuler(_Tp e[3], int order) +{ + Quat q; + + q.w = 0; + q.x = e[0]; + q.y = e[1]; + q.z = e[2]; + + EulerAngles ea = Eul_FromQuat(q, order); + + w = ea.w; + x = ea.x; + y = ea.y; + z = ea.z; +} +#endif + +template class Quaternion ; +template class Quaternion ; diff --git a/quaternion.h b/quaternion.h new file mode 100644 index 0000000..b796699 --- /dev/null +++ b/quaternion.h @@ -0,0 +1,253 @@ +//**************************************************** +//* quaternion.h * +//* * +//* Implementaion for a generalized quaternion class * +//* * +//* Written 1.25.00 by Angela Bennett * +//**************************************************** +// MODIFIED 3.14.14 by Giovanni Pinamonti + +#ifndef _QUATERNION_H_ +#define _QUATERNION_H_ + +//#include +#include +#include + +using namespace std; + +template +class Quaternion +{ + + public: + + //Quaternion + // -default constructor + // -creates a new quaternion with all parts equal to zero + Quaternion(){w=0;x=0;y=0,z=0;}; + + //Quaternion + // -constructor + // -parametes : w, x, y, z elements of the quaternion + // -creates a new quaternion based on the elements passed in + Quaternion(_Tp wi, _Tp xi, _Tp yi, _Tp zi){w=wi; x=xi; y=yi; z=zi;}; + + //Quaternion + // -constructor + // -parameters : 4D vector + // -creates a new quaternion based on the elements passed in + Quaternion(_Tp v[4]){ + w = v[0]; + x = v[1]; + y = v[2]; + z = v[3]; + }; + + //Quaternion + // -copy constructor + // -parameters : const quaternion q + // -creates a new quaternion based on the quaternion passed in + Quaternion(const Quaternion<_Tp>& q){ + w = q.w; + x = q.x; + y = q.y; + z = q.z; + }; + + //~Quaternion + // -default destructor + ~Quaternion(){}; + + //operator= + // -parameters : q1- Quaternion object + // -return values : Quaternion + // -when called on quaternion q2 sets q2 to be an object of q3 + inline Quaternion<_Tp> operator = (const Quaternion<_Tp>& q){ + w = q.w; + x = q.x; + y = q.y; + z = q.z; + return (*this); + /* questa non l'ho capita */ + }; + + //operator+ + // -parameters : q1 - Quaternion object + // -return value : Quaternion + // -when called on quaternion q2 adds q1 + q2 and returns the sum in a new quaternion + inline Quaternion<_Tp> operator + (const Quaternion<_Tp>& q){ + return Quaternion(w+q.w, x+q.x, y+q.y, z+q.z); + }; + + //operator- + // -parameters : q1- Quaternion object + // -return values : Quaternion + // -when called on q1 subtracts q1 - q2 and returns the difference as a new quaternion + inline Quaternion<_Tp> operator - (const Quaternion<_Tp>& q){ + return Quaternion(w-q.w, x-q.x, y-q.y, z-q.z); + }; + /* manca forse il -q */ + + //operator* + // -parameters : q1 - Quaternion object + // -return values : Quaternion + // -when called on a quaternion q2, multiplies q2 *q1 and returns the product in a new quaternion + inline Quaternion<_Tp> operator * (const Quaternion<_Tp>& q){ + return Quaternion( + w*q.w - x*q.x - y*q.y - z*q.z, + w*q.x + x*q.w + y*q.z - z*q.y, + w*q.y + y*q.w + z*q.x - x*q.z, + w*q.z + z*q.w + x*q.y - y*q.x); + }; + + //operator/ + // -parameters : q1 and q2- Quaternion objects + // -return values : Quaternion + // -divide q1 by q2 and returns the quotient as q1 + inline Quaternion<_Tp> operator / (Quaternion<_Tp>& q){ return ((*this) * (q.inverse()));} +; + + //operator+= + // -parameters : q1- Quaternion object + // -return values : Quaternion + // -when called on quaternion q3 adds q1 and q3 and returns the sum as q3 + inline Quaternion<_Tp>& operator += (const Quaternion<_Tp>& q){ + w += q.w; + x += q.x; + y += q.y; + z += q.z; + return (*this); +}; + + //operator-= + // -parameters : q1- Quaternion object + // -return values : Quaternion + // -when called on quaternion q3, subtracts q1 from q3 and returns the difference as q3 + inline Quaternion<_Tp>& operator -= (const Quaternion<_Tp>& q){ w += q.w; + x -= q.x; + y -= q.y; + z -= q.z; + return (*this); +}; + + //operator*= + // -parameters : q1- Quaternion object + // -return values : Quaternion + // -when called on quaternion q3, multiplies q3 by q1 and returns the product as q3 + inline Quaternion<_Tp>& operator *= (const Quaternion<_Tp>& q){ + _Tp w_val = w*q.w - x*q.x - y*q.y - z*q.z; + _Tp x_val = w*q.x + x*q.w + y*q.z - z*q.y; + _Tp y_val = w*q.y + y*q.w + z*q.x - x*q.z; + _Tp z_val = w*q.z + z*q.w + x*q.y - y*q.x; + + w = w_val; + x = x_val; + y = y_val; + z = z_val; + + return (*this); + }; + + //operator/= + // -parameters : q1- Quaternion object + // -return values : quaternion + // -when called on quaternion q3, divides q3 by q1 and returns the quotient as q3 + inline Quaternion<_Tp>& operator /= (Quaternion<_Tp>& q){ + (*this) = (*this)*q.inverse(); + return (*this); + }; + + //operator<< + // -parameters : ostream o, quaternion q + // -return values : + // -prints out a quaternion by it's components + friend inline ostream& operator << (ostream& output, const Quaternion<_Tp>& q) + { + output << "[" << q.w << ", " << "(" << q.x << ", " << q.y << ", " << q.z << ")]"; + return output; + } + + //operator!= + // -parameters : q1 and q2- Quaternion objects + // -return value : bool + // -determines if q1 and q2 and equal + inline bool operator != (const Quaternion<_Tp>& q){ + return (w!=q.w || x!=q.x || y!=q.y || z!=q.z) ? true : false; +}; + + //operator== + // -parameters : q1 and q2- Quaternion objects + // -return value : bool + // -determines if q1 and q2 and equal + inline bool operator == (const Quaternion<_Tp>& q){ + return (w==q.w && x==q.x && y==q.y && z==q.z) ? true : false; + }; + + + + + //other methods: norm, inverse, conjugate, toEuler + + //norm + // -parameters : none + // -return value : _Tp + // -when called on a quaternion object q, returns the norm of q + inline _Tp norm(){return (w*w + x*x + y*y + z*z);}; + + //magnitude + // -parameters : none + // -return value : _Tp + // -when called on a quaternion object q, returns the magnitude q + inline _Tp magnitude(){return sqrt(norm());}; + + //scale + // -parameters : s- a value to scale q1 by + // -return value: quaternion + // -returns the original quaternion with each part, w,x,y,z, multiplied by some scalar s + inline Quaternion<_Tp> scale(_Tp s){ + return Quaternion(w*s, x*s, y*s, z*s); + }; + + //inverse + // -parameters : none + // -return value : quaternion + // -when called on a quaternion object q, returns the inverse of q + inline Quaternion<_Tp> inverse(){return conjugate().scale(1/norm());}; + + //conjugate + // -parameters : none + // -return value : quaternion + // -when called on a quaternion object q, returns the conjugate of q + inline Quaternion<_Tp> conjugate(){return Quaternion(w, -x, -y, -z);}; + + //UnitQuaternion + // -parameters : none + // -return value : quaternion + // -when called on quaterion q, takes q and returns the unit quaternion of q + inline Quaternion<_Tp> UnitQuaternion(){ + return (*this).scale(1/(*this).magnitude()); + }; + + // -parameters : 3D vector of type _Tp + // -return value : void + // -when given a 3D vector, v, rotates v by the quaternion + inline void QuatRotation(_Tp v[3]){ + Quaternion <_Tp> qv(0, v[0], v[1], v[2]); + Quaternion <_Tp> qm = ( (*this) * qv ) * (*this).inverse(); + /* QUESTION: (*this) has to be normalized?? maybe */ + v[0] = qm.x; + v[1] = qm.y; + v[2] = qm.z; + }; + + + private: + // [w, (x, y, z)] + _Tp w, x, y, z; + +}; + + +#endif /* _QUATERION_H_ */ + diff --git a/trace.C b/trace.C new file mode 100644 index 0000000..1e83e23 --- /dev/null +++ b/trace.C @@ -0,0 +1,37 @@ +/* Program written by Giovanni Pinamonti + PhD student at + Scuola Internazionale Superiori di Studi Avanzati, Trieste, Italy + begin june 11th 2013 */ + +#include +#include +#include +#include +#include +#include +#include + +#include "my_malloc.h" //libraries to allocate vectors and matrices +#include "io.h" //definitions of input and output subroutines +#include "CMatrix.h" //definition of the class CMatrix + +using namespace std; + +int main(int argc, char **argv){//PROGRAM MAIN + char file_name[200]; + + if(argc>3) {cout<<"### ERROR: too many arguments ###"< 50){ + printf("Error. Reached maximum number of sweeps for Jacobi convergence. Exiting...\n"); + exit(1); + + } + + Sum=0.0; + avg_entry=0.0; + for(i=0; i < size; i++){ + for(j=i+1; j < size; j++){ + Sum += 2.0 * m[i][j]*m[i][j]; + avg_entry += 2.0*fabs(m[i][j]); + } + } + avg_entry = avg_entry/(size*(size-1.0)); + + for(p=0; p < size; p++){ + for(q=p+1; q < size; q++){ + /* For the first sweeps only perform rotations to kill the + off-diagonal matrix elements which are larger than average + */ + + if ((nsweeps <3) && (fabs(m[p][q]) < avg_entry)) continue; + + theta = (m[q][q]-m[p][p])/(2.0*m[p][q]); + + if (fabs(theta) > 1.0e+10) t = 1.0/(2*theta); + else t = 1.0/(fabs(theta) + sqrt (theta*theta+1.0)); + if (theta < 0.0) t = -t; + + c = 1.0/sqrt(t*t+1); + s = t*c; + tau = s/(1.0+c); + + /* Apply the Jacobi rotation to get the new M matrix */ + for(r=0; r < size; r++){ + rowp[r]=m[r][p]; + rowq[r]=m[r][q]; + } + + m[p][p] = rowp[p]- t*rowp[q]; + m[q][q] = rowq[q]+ t*rowq[p]; + + for(r=0; r < size; r++){ + if (r==p) continue; + if (r==q) continue; + m[r][p]= rowp[r] - s*(rowq[r] + tau*rowp[r]); + m[p][r]= m[r][p]; + m[r][q]= rowq[r] + s*(rowp[r] - tau*rowq[r]); + m[q][r]= m[r][q]; + } + m[p][q]=0.0; + m[q][p]=0.0; + Sum = Sum - 2*rowp[q]*rowp[q]; + + + /* Combine the various rotations to get the matrix of eigenvectors + */ + + for(r=0; r < size; r++){ + rowp[r]=eigenvectors[p][r]; + rowq[r]=eigenvectors[q][r]; + } + + for(r=0; r < size; r++){ + eigenvectors[p][r] = rowp[r] - s*(rowq[r] + tau*rowp[r]); + eigenvectors[q][r] = rowq[r] + s*(rowp[r] - tau*rowq[r]); + } + + if (Sum < 1.0e-12) goto end; + nrot++; + } + + } + } + end: + + + /* the diagonal elements of the rotated matrix contain the + eigenvalues */ + + for(r=0; r < size; r++){ + eigenvalues[r] = m[r][r]; + } + + + + /* Construct the inverse matrix by doing the spectral decomposition + and eliminating the null-eigenvalue space */ + + for (j=0;j0.0000001) invm[j][i]+= 1.0/eigenvalues[l]*eigenvectors[l][j]*eigenvectors[l][i]; + } + } + } + + printf("End matrix inversion. Sorting eigenvalues... \n"); + /* Now let's sort the eigenvalues (and eigenvectors) in descending + order */ + + + for(j=1; j < size; j++){ + + i = j-1; + a = eigenvalues[j]; + for(r=0; r < size; r++) rowp[r]=eigenvectors[j][r]; + + while ((i >=0) && (eigenvalues[i] < a)){ + eigenvalues[i+1] = eigenvalues[i]; + for(r=0; r < size; r++){ + eigenvectors[i+1][r]=eigenvectors[i][r]; + } + i--; + } + eigenvalues[i+1]=a; + for(r=0; r < size; r++){ + eigenvectors[i+1][r]= rowp[r]; + } + } + + free_d1t(rowp); + free_d1t(rowq); + +} +