2 McCaskill's Algorithm -- The algorithm calculates a base paring probability matrix from the input of one sequence.
4 $Id: nrutil.h,v 1.0 2005/10/20 14:22 $;
6 Copyright (C) 2005 Yasuo Tabei <tabei@cb.k.u-tokyo.ac.jp>
8 This is free software with ABSOLUTELY NO WARRANTY.
10 This library is free software; you can redistribute it and/or
11 modify it under the terms of the GNU Lesser General Public
12 License as published by the Free Software Foundation; either
13 version 2.1 of the License, or (at your option) any later version.
15 This library is distributed in the hope that it will be useful,
16 but WITHOUT ANY WARRANTY; without even the implied warranty of
17 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 Lesser General Public License for more details.
20 You should have received a copy of the GNU Lesser General Public
21 License along with this library; if not, write to the Free Software
22 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
31 #include <cstdlib> // by katoh
38 inline const T SQR(const T a) {return a*a;}
41 inline const T MAX(const T &a, const T &b)
42 {return b > a ? (b) : (a);}
44 inline float MAX(const double &a, const float &b)
45 {return b > a ? (b) : float(a);}
47 inline float MAX(const float &a, const double &b)
48 {return b > a ? float(b) : (a);}
51 inline const T MIN(const T &a, const T &b)
52 {return b < a ? (b) : (a);}
54 inline float MIN(const double &a, const float &b)
55 {return b < a ? (b) : float(a);}
57 inline float MIN(const float &a, const double &b)
58 {return b < a ? float(b) : (a);}
61 inline const T SIGN(const T &a, const T &b)
62 {return b >= 0 ? (a >= 0 ? a : -a) : (a >= 0 ? -a : a);}
64 inline float SIGN(const float &a, const double &b)
65 {return b >= 0 ? (a >= 0 ? a : -a) : (a >= 0 ? -a : a);}
67 inline float SIGN(const double &a, const float &b)
68 {return b >= 0 ? (a >= 0 ? a : -a) : (a >= 0 ? -a : a);}
71 inline void SWAP(T &a, T &b)
72 {T dum=a; a=b; b=dum;}
74 inline void nrerror(const string error_text)
75 // Numerical Recipes standard error handler
77 cerr << "Numerical Recipes run-time error..." << endl;
78 cerr << error_text << endl;
79 cerr << "...now exiting to system..." << endl;
87 int nn; // size of array. upper index is nn-1
91 explicit NRVec(int n); // Zero-based array
92 NRVec(const T &a, int n); //initialize to constant value
93 NRVec(const T *a, int n); // Initialize to array
94 NRVec(const NRVec &rhs); // Copy constructor
95 NRVec & operator=(const NRVec &rhs); //assignment
96 NRVec & operator=(const T &a); //assign a to every element
97 inline T & operator[](const int i); //i¡Çth element
98 inline const T & operator[](const int i) const;
99 void Allocator(int i);
100 inline int size() const;
105 NRVec<T>::NRVec() : nn(0), v(0) {}
108 NRVec<T>::NRVec(int n) : nn(n), v(new T[n]) {}
111 NRVec<T>::NRVec(const T& a, int n) : nn(n), v(new T[n])
113 for(int i=0; i<n; i++)
118 NRVec<T>::NRVec(const T *a, int n) : nn(n), v(new T[n])
120 for(int i=0; i<n; i++)
125 void NRVec<T>::Allocator(int n = 0)
131 NRVec<T>::NRVec(const NRVec<T> &rhs) : nn(rhs.nn), v(new T[nn])
133 for(int i=0; i<nn; i++)
138 NRVec<T> & NRVec<T>::operator=(const NRVec<T> &rhs)
139 // postcondition: normal assignment via copying has been performed;
140 // if vector and rhs were different sizes, vector
141 // has been resized to match the size of rhs
146 if (v != 0) delete [] (v);
150 for (int i=0; i<nn; i++)
157 NRVec<T> & NRVec<T>::operator=(const T &a) //assign a to every element
159 for (int i=0; i<nn; i++)
165 inline T & NRVec<T>::operator[](const int i) //subscripting
171 inline const T & NRVec<T>::operator[](const int i) const //subscripting
177 inline int NRVec<T>::size() const
197 NRMat(int n, int m); // Zero-based array
198 NRMat(const T &a, int n, int m); //Initialize to constant
199 NRMat(const T *a, int n, int m); // Initialize to array
200 NRMat(const NRMat &rhs); // Copy constructor
201 void Allocator(int n, int m);
202 void Allocator(const T &a, int n, int m);
203 void Allocator(const T *a, int n, int m);
204 NRMat & operator=(const NRMat &rhs); //assignment
205 NRMat & operator=(const T &a); //assign a to every element
206 inline T* operator[](const int i); //subscripting: pointer to row i
207 inline const T* operator[](const int i) const;
208 inline T & ref(const int i, const int j);
209 inline const T ref(const int i, const int j) const;
210 inline int nrows() const;
211 inline int ncols() const;
216 NRMat<T>::NRMat() : nn(0), mm(0), v(0) {}
219 NRMat<T>::NRMat(int n, int m) : nn(n), mm(m), v(new T*[n])
222 for (int i=1; i< n; i++)
227 NRMat<T>::NRMat(const T &a, int n, int m) : nn(n), mm(m), v(new T*[n])
239 NRMat<T>::NRMat(const T *a, int n, int m) : nn(n), mm(m), v(new T*[n])
251 void NRMat<T>::Allocator(int n, int m)
254 delete[] (v[0]); delete (v);
257 nn = n; mm = m; v = new T*[n];
260 for (int i=1; i< n; i++)
265 void NRMat<T>::Allocator(const T &a, int n, int m)
268 delete[] (v[0]); delete (v);
273 nn = n; mm = m; v = new T*[n];
284 void NRMat<T>::Allocator(const T *a, int n, int m)
287 delete[] (v[0]); delete (v);
292 nn = n; mm = m; v = new T*[n];
303 NRMat<T>::NRMat(const NRMat &rhs) : nn(rhs.nn), mm(rhs.mm), v(new T*[nn])
307 for (i=1; i< nn; i++)
309 for (i=0; i< nn; i++)
314 NRMat<T> & NRMat<T>::operator=(const NRMat<T> &rhs)
315 // postcondition: normal assignment via copying has been performed;
316 // if matrix and rhs were different sizes, matrix
317 // has been resized to match the size of rhs
321 if (nn != rhs.nn || mm != rhs.mm) {
331 for (i=1; i< nn; i++)
333 for (i=0; i< nn; i++)
341 NRMat<T> & NRMat<T>::operator=(const T &a) //assign a to every element
343 for (int i=0; i< nn; i++)
344 for (int j=0; j<mm; j++)
350 inline T* NRMat<T>::operator[](const int i) //subscripting: pointer to row i
356 inline const T* NRMat<T>::operator[](const int i) const
362 inline T & NRMat<T>::ref(const int i, const int j)
368 inline const T NRMat<T>::ref(const int i, const int j) const
374 inline int NRMat<T>::nrows() const
380 inline int NRMat<T>::ncols() const
403 NRMat3d(int n, int m, int k);
404 inline void Allocator(int n, int m, int k);
405 inline T** operator[](const int i); //subscripting: pointer to row i
406 inline const T* const * operator[](const int i) const;
407 inline int dim1() const;
408 inline int dim2() const;
409 inline int dim3() const;
414 NRMat3d<T>::NRMat3d(): nn(0), mm(0), kk(0), v(0) {}
416 NRMat3d<T>::NRMat3d(int n, int m, int k) : nn(n), mm(m), kk(k), v(new T**[n])
420 v[0][0] = new T[n*m*k];
422 v[0][j] = v[0][j-1] + k;
425 v[i][0] = v[i-1][0] + m*k;
427 v[i][j] = v[i][j-1] + k;
432 inline void NRMat3d<T>::Allocator(int n, int m, int k)
436 v[0][0] = new T[n*m*k];
438 v[0][j] = v[0][j-1] + k;
441 v[i][0] = v[i-1][0] + m*k;
443 v[i][j] = v[i][j-1] + k;
448 inline T** NRMat3d<T>::operator[](const int i) //subscripting: pointer to row i
454 inline const T* const * NRMat3d<T>::operator[](const int i) const
460 inline int NRMat3d<T>::dim1() const
466 inline int NRMat3d<T>::dim2() const
472 inline int NRMat3d<T>::dim3() const
478 NRMat3d<T>::~NRMat3d()
487 //The next 3 classes are used in artihmetic coding, Huffman coding, and
488 //wavelet transforms respectively. This is as good a place as any to put them!
491 NRVec<unsigned long> *ilob_p,*iupb_p,*ncumfq_p;
493 NRVec<unsigned long> &ilob,&iupb,&ncumfq;
494 unsigned long jdif,nc,minint,nch,ncum,nrad;
495 arithcode(unsigned long n1, unsigned long n2, unsigned long n3)
496 : ilob_p(new NRVec<unsigned long>(n1)),
497 iupb_p(new NRVec<unsigned long>(n2)),
498 ncumfq_p(new NRVec<unsigned long>(n3)),
499 ilob(*ilob_p),iupb(*iupb_p),ncumfq(*ncumfq_p) {}
501 if (ilob_p != 0) delete ilob_p;
502 if (iupb_p != 0) delete iupb_p;
503 if (ncumfq_p != 0) delete ncumfq_p;
509 NRVec<unsigned long> *icod_p,*ncod_p,*left_p,*right_p;
511 NRVec<unsigned long> &icod,&ncod,&left,&right;
513 huffcode(unsigned long n1, unsigned long n2, unsigned long n3,
515 icod_p(new NRVec<unsigned long>(n1)),
516 ncod_p(new NRVec<unsigned long>(n2)),
517 left_p(new NRVec<unsigned long>(n3)),
518 right_p(new NRVec<unsigned long>(n4)),
519 icod(*icod_p),ncod(*ncod_p),left(*left_p),right(*right_p) {}
521 if (icod_p != 0) delete icod_p;
522 if (ncod_p != 0) delete ncod_p;
523 if (left_p != 0) delete left_p;
524 if (right_p != 0) delete right_p;
530 NRVec<DP> *cc_p,*cr_p;
534 wavefilt() : cc(*cc_p),cr(*cr_p) {}
535 wavefilt(const DP *a, const int n) : //initialize to array
536 cc_p(new NRVec<DP>(n)),cr_p(new NRVec<DP>(n)),
537 ncof(n),ioff(-(n >> 1)),joff(-(n >> 1)),cc(*cc_p),cr(*cr_p) {
542 for (i=0; i<n; i++) {
548 if (cc_p != 0) delete cc_p;
549 if (cr_p != 0) delete cr_p;
554 /* Triangle Matrix Class
555 ---------------------------------------------------------
556 |v[0][0]|v[0][1]| | | |v[0][n-1]|
557 |-------|-------|------|------------|---------|---------|
558 |v[1][1]|v[1][2]| | |v[1][n-2]| |
559 |-------|-------|------|------------|---------|---------|
561 |-------|-------|------|------------|---------|---------|
565 |-------|-----------------------------------------------|
567 |-------|-----------------------------------------------|
568 |v[n-2][0]|v[n-2][1]| |
569 |-------|-----------------------------------------------|
571 |-------------------------------------------------------|
578 inline T* operator[](const int i); //subscripting: pointer to row i
579 inline const T* operator[](const int i) const;
582 Trimat(int n); // Zero-based array
583 Trimat(const T &a, int n); //Initialize to constant
584 Trimat(const T *a, int n); // Initialize to array
585 Trimat(const Trimat &rhs); // Copy constructor
586 void Allocator(int n);
587 void Allocator(const T &a, int n);
588 void Allocator(const T *a, int n);
589 Trimat & operator=(const Trimat &rhs); //assignment
590 Trimat & operator=(const T &a); //assign a to every element
591 inline T & ref(const int i, const int j);
592 inline T * getPointer(const int i, const int j);
593 inline T * begin() const;
594 inline T * end() const;
595 inline const T ref(const int i, const int j) const;
596 inline int nrows() const;
601 Trimat<T>::Trimat() : nn(0), v(0) {}
604 Trimat<T>::Trimat(int n) : nn(n), v(new T*[n])
606 v[0] = new T[n*(n+1)/2];
607 for (int i=1; i< n; i++)
608 v[i] = v[i-1] + (n-i+1);
610 for (int i=0; i< n; i++)
611 for (int j=0; j<(n-i); j++)
615 Trimat<T>::Trimat(const T &a, int n) : nn(n), v(new T*[n])
618 v[0] = new T[n*(n+1)/2];
620 v[i] = v[i-1] + (n-i+1);
622 for (j=0; j<(n-i); j++)
627 Trimat<T>::Trimat(const T *a, int n) : nn(n), v(new T*[n])
630 v[0] = new T[n*(n+1)/2];
632 v[i] = v[i-1] + (n-i+1);
634 for (j=0; j<(n-i); j++)
640 void Trimat<T>::Allocator(int n)
642 nn = n; v = new T*[n];
644 v[0] = new T[n*(n+1)/2];
645 for (int i=1; i< n; i++)
646 v[i] = v[i-1] + (n-i+1);
650 void Trimat<T>::Allocator(const T &a, int n)
652 nn = n; v = new T*[n];
655 v[0] = new T[n*(n+1)/2];
656 for (i=1; i < n; i++)
657 v[i] = v[i-1] + (n-i+1);
658 for (i=0; i < n; i++)
659 for (j=0; j < (n-i); j++)
664 void Trimat<T>::Allocator(const T *a, int n)
666 nn = n; v = new T*[n];
668 v[0] = new T[n*(n+1)/2];
670 v[i] = v[i-1] + (n-i+1);
672 for (j=0; j<(n-i); j++)
678 Trimat<T>::Trimat(const Trimat &rhs) : nn(rhs.nn), v(new T*[nn])
681 v[0] = new T[nn*(nn+1)/2];
682 for (i=1; i< nn; i++)
683 v[i] = v[i-1] + (nn-i+1);
684 for (i=0; i< nn; i++)
685 for (j=0; j<(nn-i); j++)
689 Trimat<T> & Trimat<T>::operator=(const Trimat<T> &rhs)
690 // postcondition: normal assignment via copying has been performed;
691 // if matrix and rhs were different sizes, matrix
692 // has been resized to match the size of rhs
703 v[0] = new T[nn*(nn+1)/2];
705 for (i=1; i< nn; i++)
706 v[i] = v[i-1] + (nn-i+1);
707 for (i=0; i< nn; i++)
708 for (j=0; j<(nn-i); j++)
715 Trimat<T> & Trimat<T>::operator=(const T &a) //assign a to every element
717 for (int i=0; i< nn; i++)
718 for (int j=0; j<nn-i; j++)
724 inline T & Trimat<T>::ref(const int i, const int j)
730 inline const T Trimat<T>::ref(const int i, const int j) const
736 inline T * Trimat<T>::getPointer(const int i, const int j)
742 inline T * Trimat<T>::begin() const
748 inline T * Trimat<T>::end() const
750 return (&v[nn-1][0] + 1);
754 inline int Trimat<T>::nrows() const
769 /* Triangle Vertical Matrix Class
770 ---------------------------------------------------------
771 |v[0][0]|v[1][0]| | | |v[n-1][0]|
772 |-------|-------|------|------------|---------|---------|
773 |v[0][1]|v[1][1]| | |v[n-2][1]| |
774 |-------|-------|------|------------|---------|---------|
776 |-------|-------|------|------------|---------|---------|
780 |-------|-----------------------------------------------|
782 |-------|-----------------------------------------------|
783 |v[0][n-2]|v[n-2][n-2]| |
784 |-------|-----------------------------------------------|
786 |-------------------------------------------------------|
793 inline T* operator[](const int i); //subscripting: pointer to row i
794 inline const T* operator[](const int i) const;
797 TriVertMat(int n); // Zero-based array
798 TriVertMat(const T &a, int n); //Initialize to constant
799 TriVertMat(const T *a, int n); // Initialize to array
800 TriVertMat(const TriVertMat &rhs); // Copy constructor
801 void Allocator(int n);
802 void Allocator(const T &a, int n);
803 void Allocator(const T *a, int n);
804 TriVertMat & operator=(const TriVertMat &rhs); //assignment
805 TriVertMat & operator=(const T &a); //assign a to every element
806 inline T & ref(const int i, const int j);
807 inline T * getPointer(const int i, const int j);
808 inline const T ref(const int i, const int j) const;
809 inline int nrows() const;
814 TriVertMat<T>::TriVertMat() : nn(0), v(0) {}
817 TriVertMat<T>::TriVertMat(int n) : nn(n), v(new T*[n])
819 v[0] = new T[n*(n+1)/2];
820 for (int i=1; i< n; i++)
825 TriVertMat<T>::TriVertMat(const T &a, int n) : nn(n), v(new T*[n])
828 v[0] = new T[n*(n+1)/2];
832 for (j=0; j<(n-i); j++)
837 TriVertMat<T>::TriVertMat(const T *a, int n) : nn(n), v(new T*[n])
840 v[0] = new T[n*(n+1)/2];
844 for (j=0; j<(n-i); j++)
850 void TriVertMat<T>::Allocator(int n)
852 nn = n; v = new T*[n];
854 v[0] = new T[n*(n+1)/2];
855 for (int i=1; i< n; i++)
860 void TriVertMat<T>::Allocator(const T &a, int n)
862 nn = n; v = new T*[n];
865 v[0] = new T[n*(n+1)/2];
869 for (j=0; j<(n-i); j++)
874 void TriVertMat<T>::Allocator(const T *a, int n)
876 nn = n; v = new T*[n];
878 v[0] = new T[n*(n+1)/2];
882 for (j=0; j<(n-i); j++)
888 TriVertMat<T>::TriVertMat(const TriVertMat &rhs) : nn(rhs.nn), v(new T*[nn])
891 v[0] = new T[nn*(nn+1)/2];
892 for (i=1; i< nn; i++)
894 for (i=0; i< nn; i++)
895 for (j=0; j<(nn-i); j++)
899 TriVertMat<T> & TriVertMat<T>::operator=(const TriVertMat<T> &rhs)
900 // postcondition: normal assignment via copying has been performed;
901 // if matrix and rhs were different sizes, matrix
902 // has been resized to match the size of rhs
913 v[0] = new T[nn*(nn+1)/2];
915 for (i=1; i< nn; i++)
917 for (i=0; i< nn; i++)
918 for (j=0; j<(nn-i); j++)
925 TriVertMat<T> & TriVertMat<T>::operator=(const T &a) //assign a to every element
927 for (int i=0; i< nn; i++)
928 for (int j=0; j<nn-i; j++)
934 inline T & TriVertMat<T>::ref(const int i, const int j)
940 inline const T TriVertMat<T>::ref(const int i, const int j) const
946 inline T * TriVertMat<T>::getPointer(const int i, const int j)
952 inline int TriVertMat<T>::nrows() const
958 TriVertMat<T>::~TriVertMat()
967 //Overloaded complex operations to handle mixed float and double
968 //This takes care of e.g. 1.0/z, z complex<float>
969 inline const complex<float> operator+(const double &a,
970 const complex<float> &b) { return float(a)+b; }
971 inline const complex<float> operator+(const complex<float> &a,
972 const double &b) { return a+float(b); }
973 inline const complex<float> operator-(const double &a,
974 const complex<float> &b) { return float(a)-b; }
975 inline const complex<float> operator-(const complex<float> &a,
976 const double &b) { return a-float(b); }
977 inline const complex<float> operator*(const double &a,
978 const complex<float> &b) { return float(a)*b; }
979 inline const complex<float> operator*(const complex<float> &a,
980 const double &b) { return a*float(b); }
981 inline const complex<float> operator/(const double &a,
982 const complex<float> &b) { return float(a)/b; }
983 inline const complex<float> operator/(const complex<float> &a,
984 const double &b) { return a/float(b); }
985 //some compilers choke on pow(float,double) in single precision. also atan2
986 inline float pow (float x, double y) {return pow(double(x),y);}
987 inline float pow (double x, float y) {return pow(x,double(y));}
988 inline float atan2 (float x, double y) {return atan2(double(x),y);}
989 inline float atan2 (double x, float y) {return atan2(x,double(y));}
991 #endif /* _NR_UTIL_H_ */