    001    package calhoun.util;
    002
    003    import java.io.BufferedWriter;
    004    import java.io.FileWriter;
    005    import java.io.IOException;
    006    import java.io.Writer;
    007
    008    import cern.colt.function.DoubleFunction;
    009    import cern.colt.matrix.DoubleMatrix2D;
    010    import cern.colt.matrix.impl.DenseDoubleMatrix2D;
    011    import cern.colt.matrix.linalg.EigenvalueDecomposition;
    012
    013    public class ColtUtil {
    014
    015            public static void exponentiate_real_matrix(DoubleMatrix2D X, DoubleMatrix2D Y, int degree)
    016            // Sets Y to equal exp(X); X is assumed to be real.
    017            // Y can be expressed as a convergent power series in X:
    018            // Y = I + X + X*X/2 + X*X*X/6 + X*X*X*X/24 + ...
    019            // the power series is truncated after the X^(degree) term.
    020            // Since X is real and each term above is real, Y is also real.
    021            // An eigenvector of X with eigenvalue lambda is an eigenvector of Y with eigenvalue exp(lambda)
    022            {
    023                    // S is a reall symmetric matrix, so it has real eigenvalues.
    024                    /*
    025                    EigenvalueDecomposition E = new EigenvalueDecomposition(X);
    026
    027                    Assert.a(degree >= 1);
    028
    029                    // A = V*D*V'
    030                    DoubleMatrix2D V = E.getV();
    031                    DoubleMatrix2D D = E.getD();
    032
    033                    System.out.println("The matrix D is: ");
    034                    System.out.println(ColtUtil.format(D));
    035                    System.out.println("The matrix V is: ");
    036                    System.out.println(ColtUtil.format(V));
    037
    038                    System.out.println("Real eigenvalues D is: ");
    039                    System.out.println(ColtUtil.format(E.getRealEigenvalues().toArray()));
    040                    System.out.println("Imaginary eigenvalues D is: ");
    041                    System.out.println(ColtUtil.format(E.getImagEigenvalues().toArray()));
    042                    */
    043                    /*
    044                    for (int i=0; i<4; i++) {
    045                            // only the diagonal entries:
    046                            D.setQuick(i,i,Math.exp(D.getQuick(i,i)));
    047                    }
    048                    */
    049                    /*
    050                    System.out.println("After exponentiating diagonal, The matrix D is: ");
    051                    System.out.println(ColtUtil.format(D));
    052
    053                    DoubleMatrix2D B = new DenseDoubleMatrix2D(4,4); B.assign(0);
    054                    V.zMult(D,B); // B=V*D;
    055
    056                    B.zMult(V,Y,1.0,0.0,false,true); // Y = 1.0*B*V' + 0.0*Y;
    057                    */
    058
    059                    //ColtMatrix Z;
    060
    061                    DenseDoubleMatrix2D EYE = new DenseDoubleMatrix2D(4,4);
    062                    EYE.assign(0.0);  for (int i=0; i<4; i++) { EYE.setQuick(i,i,1.0); }  // Intially sets Y to the identity matrix;
    063                    Y.assign(EYE);
    064                    DenseDoubleMatrix2D Z = new DenseDoubleMatrix2D(4,4);   Z.assign(EYE);
    065                    DenseDoubleMatrix2D W = new DenseDoubleMatrix2D(4,4);
    066                    for (int j=1; j<=degree; j++) {
    067                            W.assign(Z);
    068                            W.zMult(X,Z,1.0/((double) j),0.0,false,false);   // Z = (1/j)*W*X + 0*Z = X^j/j!
    069                            Z.zMult(EYE,Y,1.0,1.0,false,true);   // Y = Z*I + Y = Z+Y = Y+X^j/j!
    070                    }
    071
    072                    //Y = ((ColtMatrix) X).exponential(200);
    073
    074                    /*public DoubleMatrix2D zMult(DoubleMatrix2D B,
    075                                DoubleMatrix2D C,
    076                                double alpha,
    077                                double beta,
    078                                boolean transposeA,
    079                                boolean transposeB)
    080    Linear algebraic matrix-matrix multiplication; C = alpha * A x B + beta*C. C[i,j] = alpha*Sum(A[i,k] * B[k,j]) + beta*C[i,j], k=0..n-1.
    081    Matrix shapes: A(m x n), B(n x p), C(m x p).
    082    Note: Matrix shape conformance is checked after potential transpositions.
    083    */
    084
    085                    return;
    086            }
    087
    088
    089            public static void exponentiate_symmetric_matrix(DoubleMatrix2D X, DoubleMatrix2D Y)
    090            // Sets Y to equal exp(X); X is assumed to be real and symmetric.
    091            {
    092                    // S is a reall symmetric matrix, so it has real eigenvalues.
    093                    EigenvalueDecomposition E = new EigenvalueDecomposition(X);
    094
    095                    // A = V*D*V'
    096                    DoubleMatrix2D V = E.getV();
    097                    DoubleMatrix2D D = E.getD();
    098
    099                    for (int i=0; i<4; i++) {
    100                            // only the diagonal entries:
    101                            D.setQuick(i,i,Math.exp(D.getQuick(i,i)));
    102                    }
    103
    104                    DoubleMatrix2D B = new DenseDoubleMatrix2D(4,4); B.assign(0);
    105                    V.zMult(D,B); // B=V*D;
    106
    107                    B.zMult(V,Y,1.0,0.0,false,true); // Y = 1.0*B*V' + 0.0*Y;
    108
    109                    return;
    110            }
    111
    112            public static double dotProduct(double[] a, double[] b) {
    113                    int length = a.length;
    114                    Assert.a(a.length == b.length, length, " length array dotted with array of length ", b.length);
    115                    double ret = 0.0;
    116                    for(int i=0; i<length; ++i) {
    117                            ret += a[i]*b[i];
    118                    }
    119                    return ret;
    120            }
    121
    122            public static String format(double[] str) {
    123                    StringBuffer b = new StringBuffer();
    124                    for (double d : str) {
    125                            b.append(String.format("%e, ", d));
    126                    }
    127                    if(b.length() > 2) {
    128                            b.setLength(b.length() - 2);
    129                    }
    130                    return b.toString();
    131            }
    132
    133            public static void printToFile(double[] str, int numRows, int numCols) {
    134                    try {
    135                            Writer fout = new BufferedWriter(new FileWriter("test/working/crf_forwardPass_60kTEST.txt"));
    136
    137                            fout.write("Viterbi Map from CRF Forward Pass\n");
    138                            fout.write("\n");
    139                            fout.write("\t");
    140                            for (int r=0; r<numRows; r++)
    141                                    fout.write(String.format("%1$11d", r));
    142                            fout.write("\n\n");
    143
    144                            for (int c=0; c<numCols; c++)
    145                            {
    146                                    fout.write(c + ")\t");
    147                                    for (int r=0; r<numRows; r++)
    148                                    {
    149                                            double val = str[numRows * c + r];
    150                                            if (Double.isInfinite(val))
    151                                                    fout.write(String.format("%1$11.0f", 0.0f));
    152                                            else
    153                                                    //fout.write(String.format("%e\t", bestScores[pos][st], 2));
    154                                                    fout.write(String.format("%1$11.2f", val));
    155                                    }
    156                                    fout.write("\n");
    157                            }
    158                            fout.close();
    159
    160                    } catch (IOException e) {
    161                            // TODO Auto-generated catch block
    162                            e.printStackTrace();
    163                    }
    164            }
    165
    166            public static void print(double[] str, int numRows, int numCols) {
    167                    System.out.print("\n");
    168                    for (int r=0; r<numRows; r++)
    169                    {
    170                            for (int c=0; c<numCols; c++)
    171                            {
    172                                    double val = str[numRows * c + r];
    173                                    if (Double.isInfinite(val))
    174                                            val = 0.0;
    175                                    System.out.print(String.format("%e\t", val));
    176                            }
    177                            System.out.print("\n");
    178                    }
    179            }
    180
    181            public static void print(int[] str, int numRows, int numCols) {
    182                    System.out.print("\n");
    183                    for (int r=0; r<numRows; r++)
    184                    {
    185                            for (int c=0; c<numCols; c++)
    186                            {
    187                                    System.out.print(String.format("%d\t", str[numRows * c + r]));
    188                            }
    189                            System.out.print("\n");
    190                    }
    191            }
    192
    193        public static double norm(double ar[]) {
    194            double v = 0;
    195            for (int f = 0; f < ar.length; f++)
    196                v += ar[f]*ar[f];
    197            return Math.sqrt(v);
    198        }
    199
    200        public static String format(DoubleMatrix2D mi) {
    201                    StringBuffer b = new StringBuffer();
    202                    for (int i = 0; i < mi.rows(); ++i) {
    203                            for (int j = 0; j < mi.columns(); ++j) {
    204                                    b.append(String.format("%e, ", mi.getQuick(i, j)));
    205                            }
    206                            b.setLength(b.length() - 2);
    207                            b.append("; ");
    208                    }
    209                    return b.toString();
    210            }
    211
    212            /** Returns the row with the largest value in the given column */
    213            public static int maxInColumn(DoubleMatrix2D m, int col) {
    214                    double max = Double.NEGATIVE_INFINITY;
    215                    int bestRow = -1;
    216                    for(int j = 0; j<m.rows(); ++j) {
    217                            double current = m.getQuick(j, col);
    218                            if(current > max) {
    219                                    max = current;
    220                                    bestRow = j;
    221                            }
    222                    }
    223                    Assert.a(bestRow != -1);
    224                    return bestRow;
    225            }
    226
    227            public static int maxInColumn(double[] m, int numRows, int col) {
    228                    double max = Double.NEGATIVE_INFINITY;
    229                    int bestRow = -1;
    230                    int offset = col*numRows;
    231                    for(int j = 0; j<numRows; ++j) {
    232                            double current = m[offset + j];
    233                            if(current > max) {
    234                                    max = current;
    235                                    bestRow = j;
    236                            }
    237                    }
    238                    Assert.a(bestRow != -1);
    239                    return bestRow;
    240            }
    241
    242            public static double maxValueInColumn(double[] m, int numRows, int col) {
    243                    double max = Double.NEGATIVE_INFINITY;
    244                    int offset = col*numRows;
    245                    for(int j = 0; j<numRows; ++j) {
    246                            double current = m[offset + j];
    247                            if(current > max) {
    248                                    max = current;
    249                            }
    250                    }
    251                    return max;
    252            }
    253
    254            public static final DoubleFunction exp = new  DoubleFunction() {
    255                    public double apply(double a) {
    256                            return Math.exp(a);
    257                    }
    258            };
    259            public static final DoubleFunction ln = new DoubleFunction() {
    260                    public double apply(double a) {
    261                            return Math.log(a);
    262                    }
    263            };
    264            public static void scalarmultiply(DoubleMatrix2D S, final double t) {
    265                    DoubleFunction mult = new DoubleFunction() {
    266                            public double apply(double a) {
    267                                    return t*a;
    268                            }
    269                    };
    270                    S.assign(mult);
    271            }
    272
    273    }