View Javadoc

1   /*
2    * $Id: MatrixUtils.java,v 1.5 2010/01/11 21:22:46 pah Exp $
3    * 
4    * Created on 27 Nov 2008 by Paul Harrison (paul.harrison@manchester.ac.uk)
5    * Copyright 2008 Astrogrid. All rights reserved.
6    *
7    * This software is published under the terms of the Astrogrid 
8    * Software License, a copy of which has been included 
9    * with this distribution in the LICENSE.txt file.  
10   *
11   */
12  
13  package org.astrogrid.matrix;
14  
15  import org.apache.commons.math.special.Gamma;
16  
17  import no.uib.cipr.matrix.AGDenseMatrix;
18  import no.uib.cipr.matrix.DenseVector;
19  import no.uib.cipr.matrix.MatrixEntry;
20  import no.uib.cipr.matrix.Vector;
21  import no.uib.cipr.matrix.VectorEntry;
22  
23  /**
24   * Utility functions that operate on matrices. These are intended to operate in the same fashion as the similarly named MATLAB functions.
25   * 
26   * @author Paul Harrison (paul.harrison@manchester.ac.uk) 4 Aug 2009
27   * @version $Name:  $
28   * @since AIDA Stage 1
29   */
30  public class MatrixUtils {
31  
32      /**
33       * Forms the covariance matrix. For matrices, where each row is an
34       * observation, and each column a variable, COV(X) is the covariance matrix.
35       * DIAG(COV(X)) is a vector of variances for each column, and
36       * SQRT(DIAG(COV(X))) is a vector of standard deviations
37       * 
38       * COV(X) or COV(X,Y) normalizes by (N-1) if N>1, where N is the number of
39       * observations. This makes COV(X) the best unbiased estimate of the
40       * covariance matrix if the observations are from a normal distribution. For
41       * N=1, COV normalizes by N.
42       * 
43       * compatible with MatLab definition of the same
44       * 
45       * also @see http://mathworld.wolfram.com/Covariance.html
46       * 
47       * @param x
48       * @return
49       */
50      public static Matrix cov(Matrix x) {
51  
52          int m = x.numRows();
53          int n = x.numColumns();
54  
55          if (m == 1) {
56              return zeros(n);
57          } else {
58              // form mean of each column in y
59              Vector y = sum(x, 1);
60              y.scale(1.0 / m);
61              Matrix result = new AGDenseMatrix(m, n);
62              Matrix result2 = new AGDenseMatrix(m,n);
63              for (MatrixEntry me : x) {
64                  int col = me.column();
65                  int row = me.row();
66                  double value = me.get() - y.get(col);
67                  result.set(row, col, value); // remove the
68                                                                  // mean from
69                  result2.set(row, col, value);                                                // each value
70  
71              }
72              Matrix cov = new AGDenseMatrix(n,n);
73              result.transAmult(1.0/(m-1), result2, cov);//m-1 is the most likely use in matlab - can be overwritten there though..
74              
75              return cov;
76              
77          }
78  
79      }
80  
81      public static AGDenseMatrix zeros(int n) {
82          return (AGDenseMatrix) new AGDenseMatrix(n, n).zero();
83      }
84      public static AGDenseMatrix zeros(int n, int m) {
85          return (AGDenseMatrix) new AGDenseMatrix(n, m).zero();
86      }
87      
88      public static Matrix ones(int n) {
89          return new AGDenseMatrix(n, n).ones();
90      }
91      public static Matrix ones(int n, int m) {
92          return new AGDenseMatrix(n, m).ones();
93      }
94  
95      public static Vector sum(Matrix x, int dim) {
96          return x.sum(dim);
97      }
98      
99      public static double sum(Vector v){
100         double s = 0;
101         for (VectorEntry vectorEntry : v) {
102             s += vectorEntry.get();
103         }
104         return s;
105     }
106     
107     public static Vector sum(Matrix m) {
108         return sum(m,1);
109     }
110 
111     
112     public static Vector diag(no.uib.cipr.matrix.Matrix x){
113         if(x.isSquare()){
114             Vector retval = new DenseVector(x.numColumns());
115             for (int i = 0; i < retval.size(); i++) {
116                 retval.set(i, x.get(i, i));
117             }
118             return retval;
119         } else {
120             throw new IllegalArgumentException("the matrix must be square");
121         }
122     }
123     
124     public static Matrix diag(Vector v){
125         Matrix retval = new AGDenseMatrix(v.size(), v.size());
126         for (int i = 0; i < v.size(); i++) {
127             retval.set(i, i, v.get(i));
128         }
129         return retval;
130     }
131     
132     public static Matrix repmatv(Vector v, int ni, int nj)
133     {
134         Matrix retval = new AGDenseMatrix(v, ni, nj);
135         return retval;
136     }
137     public static Matrix repmat(Matrix v, int ni, int nj)
138     {
139         Matrix retval = new AGDenseMatrix(v, ni, nj);
140         return retval;
141     }
142     
143     /**
144      * Repeat a vector into a matrix treating the vector as a column vector.
145      * @param v
146      * @param ni
147      * @param nj
148      * @return
149      */
150      public static Matrix repmat(Vector v, int ni, int nj)
151     {
152         Matrix retval = new AGDenseMatrix(v);
153         return repmat(retval, ni, nj);
154         
155     }
156     
157     /**
158      * Repeat a vector into a matrix treating the vector as a row vector.
159      * @param v
160      * @param ni
161      * @param nj
162      * @return
163      */
164     public static Matrix repmatt(Vector v, int ni, int nj)
165     {
166         Matrix retval = transpose(new AGDenseMatrix(v));
167         return repmat(retval, ni, nj);
168         
169     }
170      
171     public static AGDenseMatrix reshape(Matrix m, int ni, int nj){
172         AGDenseMatrix retval = new AGDenseMatrix(m);
173         return retval.reshape(ni, nj);
174     }
175     
176     public static AGDenseMatrix reshape(Vector v, int ni, int nj ){
177         AGDenseMatrix retval = new AGDenseMatrix(v);
178         return  retval.reshape(ni, nj);
179     }
180     
181     public static double det(Matrix m){
182         return m.det();
183     }
184     
185     public static double trace(Matrix m){
186         return m.trace();
187     }
188     
189     /**
190      * Elementwise divide of matrix.
191      * @param a
192      * @param b
193      * @return
194      */
195     public static Matrix divide(Matrix  a, Matrix b){
196        
197          if(a.numRows() == b.numRows() && a.numColumns() == b.numColumns()){ 
198          Matrix retval = new AGDenseMatrix(a.numRows(), a.numColumns());
199          a.divide(b, retval);
200          return retval;
201          }
202          else {
203              throw new IllegalArgumentException("matrices are not same size");
204          }
205     }
206     /**
207      * Elementwise multiply of matrix.
208      * @param a
209      * @param b
210      * @return
211      */
212     public static Matrix times(no.uib.cipr.matrix.Matrix  a, no.uib.cipr.matrix.Matrix b){
213        
214          if(a.numRows() == b.numRows() && a.numColumns() == b.numColumns()){ 
215          Matrix retval = new AGDenseMatrix(a.numRows(), a.numColumns());
216          for (int i = 0; i < retval.numRows(); i++) {
217             for (int j = 0; j < retval.numColumns(); j++) {
218                 retval.set(i, j, a.get(i, j)*b.get(i, j));
219             }
220         }
221          return retval;
222          }
223          else {
224              throw new IllegalArgumentException("matrices are not same size");
225          }
226     }
227     /**
228      * Elementwise divide of vector.
229      * @param a
230      * @param b
231      * @return
232      */
233     public static Vector divide(Vector  a, Vector b){
234        
235          if(a.size() == b.size() ){ 
236          DenseVector retval = new DenseVector(a.size());
237          for (int i = 0; i < retval.size(); i++) {
238             retval.set(i, a.get(i)/b.get(i));
239         }
240          return retval;
241          }
242          else {
243              throw new IllegalArgumentException("vectors are not same size");
244          }
245     }
246     /**
247      * Elementwise addition of vectors.
248      * @param a
249      * @param b
250      * @return
251      */
252     public static Vector add(Vector  a, Vector b){
253        
254          if(a.size() == b.size() ){ 
255          DenseVector retval = new DenseVector(a.size());
256          for (int i = 0; i < retval.size(); i++) {
257             retval.set(i, a.get(i)+b.get(i));
258         }
259          return retval;
260          }
261          else {
262              throw new IllegalArgumentException("vectors are not same size");
263          }
264     }
265     /**
266      * Elementwise product of vectors.
267      * @param a
268      * @param b
269      * @return
270      */
271     public static Vector times(Vector  a, Vector b){
272        
273          if(a.size() == b.size() ){ 
274          DenseVector retval = new DenseVector(a.size());
275          for (int i = 0; i < retval.size(); i++) {
276             retval.set(i, a.get(i)*b.get(i));
277         }
278          return retval;
279          }
280          else {
281              throw new IllegalArgumentException("vectors are not same size");
282          }
283     }
284     /**
285      * Elementwise product of vectors and double
286      * @param a
287      * @param b
288      * @return
289      */
290     public static Vector times(Vector  a, double b){
291        
292          
293          DenseVector retval = new DenseVector(a);
294          
295          return retval.scale(b);
296      }
297     /**
298      * Elementwise product of vectors and double
299      * @param a
300      * @param b
301      * @return
302      */
303     public static Vector times( double b, Vector  a){
304        
305          
306          DenseVector retval = new DenseVector(a);
307          
308          return retval.scale(b);
309      }
310     
311     public static Matrix times(double b, Matrix a){
312         AGDenseMatrix retval = new AGDenseMatrix(a);
313         return (Matrix) retval.scale(b);
314     }
315     public static Matrix times(Matrix a,double b){
316         AGDenseMatrix retval = new AGDenseMatrix(a);
317         return (Matrix) retval.scale(b);
318     }
319    
320    /**
321      * Elementwise difference of vectors.
322      * @param a
323      * @param b
324      * @return
325      */
326     public static Vector sub(Vector  a, Vector b){
327        
328          if(a.size() == b.size() ){ 
329          DenseVector retval = new DenseVector(a.size());
330          for (int i = 0; i < retval.size(); i++) {
331             retval.set(i, a.get(i)-b.get(i));
332         }
333          return retval;
334          }
335          else {
336              throw new IllegalArgumentException("vectors are not same size");
337          }
338     }
339     
340     public static Matrix inv(Matrix a){
341         return ((Matrix)a).inv();
342     }
343     
344     /**
345      * C = AB
346      * @param a
347      * @param b
348      * @return
349      */
350     public static AGDenseMatrix mult(no.uib.cipr.matrix.Matrix a, no.uib.cipr.matrix.Matrix b){
351         AGDenseMatrix c = new AGDenseMatrix(a.numRows(), b.numColumns());
352         a.mult(b, c);
353         return c;
354     }
355     /**
356      * C = AB
357      * @param a
358      * @param b
359      * @return
360      */
361     public static AGDenseMatrix mult(no.uib.cipr.matrix.Matrix a, Vector b){
362         AGDenseMatrix c = new AGDenseMatrix(a.numRows(), 1);
363         a.mult(new AGDenseMatrix(b), c);
364         return c;
365     }
366     
367   /**
368      * C = AB<sup>T</sup>
369      * @param a
370      * @param b
371      * @return
372      */
373     public static AGDenseMatrix multBt(no.uib.cipr.matrix.Matrix a, no.uib.cipr.matrix.Matrix b) {
374         Matrix c = new AGDenseMatrix(a.numRows(), b.numRows());
375         
376         return (AGDenseMatrix) a.transBmult(b, c);
377     }
378     
379     /**
380      * C = A<sup>T</sup>B
381      * @param a
382      * @param b
383      * @return
384      */
385     public static AGDenseMatrix multAt(no.uib.cipr.matrix.Matrix a, no.uib.cipr.matrix.Matrix b) {
386         Matrix c = new AGDenseMatrix(a.numColumns(), b.numColumns());
387         
388         return (AGDenseMatrix) a.transAmult(b, c);
389     }
390     
391     
392     /**
393      * C = A*B*A<sup>T</sup>
394      * @param a
395      * @param b
396      * @return
397      */
398     public static Matrix multABAT(no.uib.cipr.matrix.Matrix a, no.uib.cipr.matrix.Matrix b)
399     {
400         Matrix c = new AGDenseMatrix(a.numRows(), b.numColumns());
401         a.mult(b, c);
402         Matrix d = new AGDenseMatrix(a.numRows(), a.numRows());
403         c.transBmult(a, d);
404         return d;
405     }
406     /**
407      * C = A<sup>T</sup>*B*A
408      * @param a
409      * @param b
410      * @return
411      */
412     public static Matrix multATBA(DenseVector v, no.uib.cipr.matrix.Matrix b)
413     {
414         Matrix a = new AGDenseMatrix(new double[][]{v.getData()}); //this transposes the vector to be a row vector
415         Matrix c = new AGDenseMatrix(1, b.numColumns());
416         a.mult(b, c);
417         Matrix d = new AGDenseMatrix(1,1);
418         c.transBmult(a, d);
419         return d;
420     }
421     /**
422      * a*b - treating a as a column vector.
423      * @param a
424      * @param b
425      * @return
426      */
427     public static Matrix mult(Vector a, no.uib.cipr.matrix.Matrix  b){
428         Matrix c = new AGDenseMatrix(a);
429         Matrix d = new AGDenseMatrix(a.size(),b.numColumns());
430         return (Matrix) c.mult(b, d);
431     }
432     /**
433      * C = A<sup>T</sup>B. Make the vector  row vector before multiplying.
434      * @param a
435      * @param b
436      * @return
437      */
438     public static AGDenseMatrix multAt(Vector v, no.uib.cipr.matrix.Matrix b) {
439         Matrix a = new AGDenseMatrix(v);
440         Matrix c = new AGDenseMatrix(a.numColumns(), b.numColumns());
441         return (AGDenseMatrix) a.transAmult(b, c);
442     }
443   
444     /**
445      * R = A+B. Elementwise addition of matrices - must be the same size.
446      * @param a
447      * @param b
448      * @return
449      */
450     public static AGDenseMatrix add(no.uib.cipr.matrix.Matrix a, no.uib.cipr.matrix.Matrix b){
451         AGDenseMatrix c = new AGDenseMatrix(a);
452         return   (AGDenseMatrix) c.add(b);
453         
454     }
455     
456     public static Vector add(Vector v, double d){
457         Vector retval = new DenseVector(v.size());
458         for (int i = 0; i < v.size(); i++) {
459             retval.set(i, v.get(i)+d);
460         }
461         return retval;
462     }
463     
464     public static Vector add(double d, Vector v ){
465         return add(v,d);
466     } 
467     
468     public static Matrix add(Matrix q, double eps) {
469         AGDenseMatrix retval = new AGDenseMatrix(q.numRows(), q.numColumns());
470         for (MatrixEntry matrixEntry : q) {
471             retval.set(matrixEntry.row(), matrixEntry.column(), matrixEntry.get()+eps);
472         }
473         return retval;
474     }
475     public static Matrix add(double eps, Matrix q) {
476         return add(q,eps);
477     }
478 
479     
480     /**
481      * R = A-B. Elementwise subtraction of matrices - must be the same size.
482      * @param a
483      * @param b
484      * @return
485      */
486    public static Matrix sub(no.uib.cipr.matrix.Matrix a, no.uib.cipr.matrix.Matrix b){
487         Matrix c = new AGDenseMatrix(a);
488         return (Matrix) c.add(-1.0, b);
489         
490     }
491    public static Matrix sub(no.uib.cipr.matrix.Matrix a, double b){
492        Matrix c = new AGDenseMatrix(a);
493        
494        for (MatrixEntry matrixEntry : c) {
495         matrixEntry.set(matrixEntry.get()-b);
496     }
497        return c;
498        
499    }
500    
501    /**
502     * R= d-v = elementwise subtraction.
503  * @param d
504  * @param v
505  * @return
506  */
507 public static Vector sub(double d, Vector v){
508       Vector retval = new DenseVector(v.size());
509        for (int i = 0; i < v.size(); i++) {
510            retval.set(i, d-v.get(i));
511        }
512        return retval;
513    }
514 public static Vector sub( Vector v,double d){
515     Vector retval = new DenseVector(v.size());
516      for (int i = 0; i < v.size(); i++) {
517          retval.set(i, v.get(i)-d);
518      }
519      return retval;
520  }
521     
522     /**
523      * produces square diagonal matrix using val.
524      * @param ndim
525      * @param val - the value to put on the diagonal.
526      * @return
527      */
528     public static Matrix eye(int ndim, double val){
529        Matrix retval = new AGDenseMatrix(ndim,ndim);
530        retval.zero();
531        for (int i = 0; i < ndim; i++) {
532         retval.set(i, i, val);
533     }
534        return retval;
535     }
536     /**
537      * produces square identity matrix.
538      * @param ndim - the dimension of the matrix
539      * @return
540      */
541     public static Matrix eye(int ndim){
542        return eye(ndim,1.0);
543      }
544    
545     
546     public static AGDenseMatrix transpose(no.uib.cipr.matrix.Matrix a){
547         Matrix retval = new AGDenseMatrix(a.numColumns(), a.numRows());
548         return  (AGDenseMatrix) a.transpose(retval);
549     }
550     
551     public static double max(Vector iv) {
552        double retval = iv.get(0);
553        for (int i = 1; i < iv.size(); i++) {
554         if (iv.get(i) > retval) {
555             retval = iv.get(i);
556         }
557     }
558        return retval;
559     }
560  
561     public static Vector max(Matrix im) {
562         Vector retval = new DenseVector(im.numColumns());
563         for (int i = 0; i < im.numColumns(); i++) {
564             double maxval=im.get(0, i);
565             for (int j = 1; j < im.numRows(); j++) {
566                 if (im.get(j,i) > maxval) {
567                     maxval = im.get(j,i);
568                 }
569         
570             }
571             retval.set(i, maxval);
572      }
573         return retval;
574      }
575 
576     
577     public static Matrix rand (int i, int j) {
578         Matrix retval = new AGDenseMatrix(i, j);
579         for (MatrixEntry matrixEntry : retval) {
580             matrixEntry.set(Math.random());
581         }
582         return retval;
583     }
584     public static Matrix dirichlet_sample (Vector m, int j) {
585         Matrix retval = new AGDenseMatrix(m.size(), j);
586         for (MatrixEntry matrixEntry : retval) {
587             matrixEntry.set(Math.random()); //FIXME - this is not acutually returning dirichlet sampled random numbers.
588         }
589         return retval;
590     }
591     
592     
593     public static Matrix psi(Matrix m) {
594         
595         for (MatrixEntry matrixEntry : m) {
596             matrixEntry.set(Gamma.digamma(matrixEntry.get()));
597         }
598         return m;
599     }
600     
601     public static double psi(double d){
602         return Gamma.digamma(d);
603     }
604     
605     public static Matrix log(Matrix m) {
606         Matrix retval = new AGDenseMatrix(m.numRows(), m.numColumns());
607         for (MatrixEntry matrixEntry : m) {
608             retval.set(matrixEntry.row(),matrixEntry.column(), Math.log(matrixEntry.get()));
609         }
610         return retval;
611     }
612     public static Matrix exp(Matrix m) {
613         Matrix retval = new AGDenseMatrix(m.numRows(), m.numColumns());
614         for (MatrixEntry matrixEntry : m) {
615             retval.set(matrixEntry.row(),matrixEntry.column(), Math.exp(matrixEntry.get()));
616         }
617         return retval;
618     }
619     
620     public static Vector exp(Vector v){
621         Vector retval = new DenseVector(v.size());
622         for (int i = 0; i < v.size(); i++) {
623             retval.set(i, Math.exp(v.get(i)));
624         }
625         return retval;
626         
627     }  public static Vector log(Vector v){
628         Vector retval = new DenseVector(v.size());
629         for (int i = 0; i < v.size(); i++) {
630             retval.set(i, Math.log(v.get(i)));
631         }
632         return retval;
633         
634     }
635     /**
636      * raise each member to the power.
637      * @param m
638      * @param exp
639      * @return
640      */
641     public static Matrix pow(Matrix m, double exp){
642         Matrix retval = new AGDenseMatrix(m.numRows(), m.numColumns());
643         for (MatrixEntry matrixEntry : m) {
644             retval.set(matrixEntry.row(),matrixEntry.column(), Math.pow(matrixEntry.get(), exp));
645         }
646         return retval;
647     }
648     public static Vector pow(Vector v, double exp) {
649         Vector retval = new DenseVector(v.size());
650         for (int i = 0; i < v.size(); i++) {
651             retval.set(i, Math.pow(v.get(i),exp));
652         }
653         return retval;
654     }
655 
656     
657     public static Matrix recip(Matrix m){
658         Matrix retval = new AGDenseMatrix(m.numRows(), m.numColumns());
659         for (MatrixEntry matrixEntry : m) {
660             retval.set(matrixEntry.row(),matrixEntry.column(), 1.0/matrixEntry.get());
661         }
662         return retval;
663        
664     }
665     public static Vector recip(Vector v) {
666         Vector retval = new DenseVector(v.size());
667         for (int i = 0; i < v.size(); i++) {
668             retval.set(i, 1.0/v.get(i));
669         }
670         return retval;
671     }
672 
673     /**
674      * Form the elementwise a/m<sub>ij</sub>.
675      * @param m
676      * @param a
677      * @return
678      */
679     public static Matrix recip( double a, Matrix m){
680         Matrix retval = new AGDenseMatrix(m.numRows(), m.numColumns());
681         for (MatrixEntry matrixEntry : m) {
682             retval.set(matrixEntry.row(),matrixEntry.column(), a/matrixEntry.get());
683         }
684         return retval;
685        
686     }
687     /**
688      * Form the elementwise a/v<sub>i</sub>.
689      * @param m
690      * @param a
691      * @return
692      */
693    public static Vector recip(double a, Vector v) {
694         Vector retval = new DenseVector(v.size());
695         for (int i = 0; i < v.size(); i++) {
696             retval.set(i, a/v.get(i));
697         }
698         return retval;
699     }
700 
701     
702     /**
703      * produces a matrix which is the sum of the squares of each row.
704      * @param m
705      * @return
706      */
707     public static Matrix sumsq(Matrix m)
708     {
709         Matrix retval = new AGDenseMatrix(m.numRows(),1);
710         for (MatrixEntry matrixEntry : m) {
711              Double val = matrixEntry.get();
712              int row = matrixEntry.row();
713              retval.set(row, 1, retval.get(row, 1) + val* val);
714         }
715         return retval;
716     }
717     
718     
719     /**
720      * Return a vector filled with a sequence from 1 to n;
721      * @param n
722      * @return
723      */
724     public static Vector seq(int n){
725         Vector retval = new DenseVector(n);
726         int j = 1;
727         for (VectorEntry vectorEntry : retval) {
728             vectorEntry.set(j++);
729         }
730         return retval;
731     }
732     
733     
734     /**
735      * Forms the vector product of the two vectors -i.e. ab<sup>T</sup>.
736      * @param a
737      * @param b
738      * @return
739      */
740     public static Matrix vprod(Vector a, Vector b){
741         return multBt(new AGDenseMatrix(a), new AGDenseMatrix(b));
742     }
743     
744     /**
745      * Forms the vector product of the two vectors -i.e. aa<sup>T</sup>.
746      * @param a
747      * @param b
748      * @return
749      */
750     public static Matrix vprod(Vector a){
751         return multBt(new AGDenseMatrix(a), new AGDenseMatrix(a));
752     }
753    
754     
755     public static DenseVector delElement(DenseVector v, int i){
756         if(i>= v.size() || i < 0){
757             throw new IllegalArgumentException("element to delete - out of range");
758         }
759         int size = v.size();
760         double[] newarr= new double[size-1];
761         
762         if(i != 0){
763                 // copy first half
764             System.arraycopy(v.getData(), 0, newarr, 0, i);
765         }
766         //copy second half
767         if(i < size -1){
768         System.arraycopy(v.getData(), i+1, newarr, i, size -i-1);
769         }
770         return new DenseVector(newarr, false);
771         
772     }
773     
774     public static int min(Vector v){
775         int midx = 0, idx=0;
776         double val = v.get(idx);
777         while(++idx < v.size()){
778             if(v.get(idx) < val){
779                 midx = idx;
780                 val = v.get(idx);
781             }
782         }
783         return midx;
784     }
785 
786 }
787 
788 
789 
790 
791 /*
792  * $Log: MatrixUtils.java,v $
793  * Revision 1.5  2010/01/11 21:22:46  pah
794  * reasonable numerical stability and fidelity to MATLAB results achieved
795  *
796  * Revision 1.4  2010/01/05 21:27:57  pah
797  * add delete colum/row functions
798  *
799  * Revision 1.3  2009/09/22 07:04:16  pah
800  * daily checkin
801  *
802  * Revision 1.2  2009/09/14 19:08:43  pah
803  * code runs clustering, but not giving same results as matlab exactly
804  *
805  * Revision 1.1  2009/09/07 16:06:12  pah
806  * initial transcription of the core
807  *
808  */