0001: // -- @(#)Matrix.java 1.4 99/03/22
0002:
0003: /*
0004: * This version by David G. Melvin. Originally based on Matrix.java, v1.3, by
0005: * Akio Utsugi, but extended to include linear algebra stuff.
0006: */
0007:
0008: package uk.org.ponder.matrix;
0009:
0010: import java.io.IOException;
0011: import java.io.Writer;
0012: import java.text.DecimalFormat;
0013:
0014: import uk.org.ponder.stringutil.CharWrap;
0015: import uk.org.ponder.util.Logger;
0016: import uk.org.ponder.util.UniversalRuntimeException;
0017:
0018: class MatrixStorage {
0019: double value[];
0020:
0021: int references = 1;
0022:
0023: MatrixStorage(int row, int col) {
0024: value = new double[row * col];
0025: }
0026:
0027: }
0028:
0029: /**
0030: * A class for manipulation of numerical matrix objects with basic linear
0031: * algebra and eigen system methods.
0032: *
0033: * @version @(#)Matrix.java 1.4 99/03/22, 17 Sep 1996
0034: * @author David G. Melvin
0035: */
0036:
0037: public class Matrix implements Cloneable {
0038:
0039: final static double small = 1.0e-20;
0040: final static int JacobiIter = 50;
0041:
0042: int rows, cols;
0043: MatrixStorage storage;
0044:
0045: /**
0046: * Generalised storage system: These variables allow us to map the matrix
0047: * dimensions onto the 1d array in an arbitrary manner.
0048: */
0049:
0050: int xstep, ystep;
0051: int topleft = 0;
0052:
0053: /**
0054: * Private cache of values that are relatively expensive to compute or that
0055: * are generated in tandem or likely to be used frequently after they have
0056: * been computed. NB. ensure that these are reset to null when any elements of
0057: * the matrix are changed. The private routine clearCache is provided for this
0058: * purpose. Ensure that this is always in line with the set of cached
0059: * components.
0060: */
0061:
0062: LUDecomp LUD = null;
0063: Matrix eigenVec = null;
0064: Matrix eigenVal = null;
0065:
0066: void clearCache() {
0067: LUD = null;
0068: eigenVec = null;
0069: eigenVal = null;
0070: }
0071:
0072: public static final UnaryFunction Exp = new UnaryFunction() {
0073: public double apply(double arg) {
0074: return Math.exp(arg);
0075: }
0076: };
0077: public static final UnaryFunction Sqrt = new UnaryFunction() {
0078: public double apply(double arg) {
0079: return Math.sqrt(arg);
0080: }
0081: };
0082:
0083: // The matrix has been modified, and requires its storage reallocated.
0084: // We take the opportunity to save storage by compacting rows, and
0085: // make use of System.arraycopy to blag as many elements as possible that
0086: // occur without gaps. In order to get maximum use out of System.arraycopy,
0087: // we do not guarantee that the resulting matrix is not tranposed relative
0088: // to its storage, i.e. it may be that ystep = 1 rather than xstep = 1.
0089:
0090: void reallocate() {
0091: // System.out.println("Matrix reallocated: ");
0092: // debug();
0093: double[] value = storage.value;
0094: MatrixStorage newstorage = new MatrixStorage(rows, cols);
0095: double[] newvalue = newstorage.value;
0096:
0097: // The normal case - the elements in the old matrix form contiguous rows
0098: // without gaps between them, although they may be transposed.
0099: // We may blag the entire matrix at one go.
0100:
0101: if ((xstep == 1 || ystep == 1) && ystep == cols) {
0102: // System.out.println("Case 1");
0103: System.arraycopy(value, topleft, newvalue, 0, rows * cols);
0104: }
0105:
0106: else {
0107: if (xstep == 1) {
0108: // the slightly less normal case - elements in
0109: // old matrix are the right way
0110: // up, but with gaps.
0111: // System.out.println("Case 2");
0112: for (int row = 0; row < rows; ++row) {
0113: System.arraycopy(value, topleft + row * ystep,
0114: newvalue, 0 + row * cols, cols);
0115: }
0116: ystep = cols;
0117: }
0118:
0119: else if (ystep == 1) {
0120: // the even less normal case - elements in old matrix are the wrong way
0121: // up, and with gaps.
0122: // System.out.println("Case 3");
0123: for (int col = 0; col < cols; ++col) {
0124: System.arraycopy(value, topleft + col * xstep,
0125: newvalue, 0 + col * rows, rows);
0126: }
0127: xstep = rows;
0128: }
0129:
0130: else {
0131: // the completely insane case - elements are arbitrarily distributed
0132: // round
0133: // the matrix with gaps everywhere. This is quite likely never to occur!
0134: int newaddress = 0;
0135: for (int row = 0; row < rows; ++row) {
0136: for (int col = 0; col < cols; ++col) {
0137: newvalue[newaddress++] = getMval(row, col);
0138: }
0139: }
0140: xstep = 1;
0141: ystep = cols;
0142: }
0143: }
0144:
0145: // We escaped alive!
0146: storage.references--;
0147: if (storage.references == 0)
0148: throw new RuntimeException(
0149: "Programming error: no references remain to MatrixStorage");
0150:
0151: storage = newstorage;
0152: topleft = 0;
0153: // System.out.println("Reallocated Matrix is now");
0154: // debug();
0155: }
0156:
0157: /**
0158: * This package private method MUST be called before any operation that would
0159: * alter the matrix contents
0160: */
0161: void aboutToModify() {
0162: clearCache();
0163: if (storage.references > 1)
0164: reallocate();
0165: }
0166:
0167: /**
0168: * This is an experimental method, noting that *ALL* of the matrix values in
0169: * storage are about to be overwritten, so we may make the optimisation of
0170: * assuring that the storage is not transposed. We may then get the benefits
0171: * of sharing storage where it is the same size, and also of assuming that
0172: * xstep = 1, and we may also step through the array from TL to BR by adding
0173: * 1, etc.
0174: */
0175: void aboutToObliterate() {
0176: clearCache();
0177: if (storage.references > 1)
0178: reallocate();
0179: if (ystep == 1) {
0180: int temp = ystep;
0181: ystep = xstep;
0182: xstep = temp;
0183: }
0184: }
0185:
0186: /** An extremely dangerous method to get the underlying array storing
0187: * the data for this matrix, for clients who have good expectations it
0188: * is in a reasonable arrangement.
0189: * @return
0190: */
0191: public double[] getStorageDirty() {
0192: return storage.value;
0193: }
0194:
0195: /**
0196: * The implementation of the matrix maps its elements into a 1d linear array.
0197: * This private interface routine computes the correct location in the 1d
0198: * array as if it were a 2d array. The advantage is the reduction in overhead
0199: * imposed by java's 2d array implementation. This also allows arbitrary
0200: * linear mapping of the matrix into 1d, allowing for transpose, submatrices,
0201: * and copies to be efficiently implemented.
0202: * <p>
0203: * For efficiency, calls to this method could be eliminated by use of a
0204: * strength-reduced index, except where this would cause the resulting code to
0205: * become unreadable.
0206: */
0207: double getMval(int y, int x) {
0208: return storage.value[topleft + x * xstep + y * ystep];
0209: }
0210:
0211: /**
0212: * The implementation of the matrix maps its elements into a 1d linear array.
0213: * This private interface routine computes the correct location in the 1d
0214: * array as if it were a 2d array. The advantage is the reduction in overhead
0215: * imposed by java's 2d array implementation. This also allows arbitrary
0216: * linear mapping of the matrix into 1d, allowing for transpose, submatrices,
0217: * and copies to be efficiently implemented.
0218: * <p>
0219: * For efficiency, calls to this method could be eliminated by use of a
0220: * strength-reduced index, except where this would cause the resulting code to
0221: * become unreadable.
0222: * <p>
0223: * Since this function is package private, it does not clear the caches for
0224: * extra speed.
0225: */
0226:
0227: void setMval(int y, int x, double val) {
0228: storage.value[topleft + x * xstep + y * ystep] = val;
0229: }
0230:
0231: /**
0232: * The column separator used when rendering a matrix into a string.
0233: */
0234: public String fieldSeparator = " ";
0235:
0236: /**
0237: * The number format to be used when rendering the matrix into a string. The
0238: * default is to emulate the behaviour of Matlab to a degree by allowing 4
0239: * digits after the decimal point.
0240: */
0241: public static DecimalFormat decFormat = new DecimalFormat(
0242: " 0.####;-0.####");
0243:
0244: /**
0245: * Construct a square matrix object, with initial values set to zero.
0246: *
0247: * @param n width and height.
0248: */
0249: public Matrix(int n) {
0250: this (n, n);
0251: }
0252:
0253: /**
0254: * Construct a rectangular matrix object, with initial values set to zero.
0255: *
0256: * @param row number of row
0257: * @param col number of column
0258: */
0259: public Matrix(int rows, int cols) {
0260: this (rows, cols, true);
0261: }
0262:
0263: /**
0264: * Private constructor which allows the choice not to allocate storage. If
0265: * storage is shared with another Matrix, remember to adjust the reference
0266: * count of the storage.
0267: *
0268: * @param row number of row
0269: * @param col number of column
0270: */
0271:
0272: Matrix(int rows, int cols, boolean allocate) {
0273: this .xstep = 1;
0274: this .ystep = cols;
0275: this .rows = rows;
0276: this .cols = cols;
0277: if (allocate)
0278: storage = new MatrixStorage(rows, cols);
0279: }
0280:
0281: /**
0282: * Construct a rectangular matrix object from a 2D array of doubles.
0283: *
0284: * @param data
0285: */
0286:
0287: public Matrix(double[][] data) {
0288: this (data, 0, 0);
0289: }
0290:
0291: /** Constructs a matrix from a multidimensional double array, with
0292: * specified space for padding.
0293: * @param data
0294: * @param padrow
0295: * @param padcol
0296: */
0297: public Matrix(double[][] data, int padrow, int padcol) {
0298: this (data.length + padrow, (data.length == 0 ? 0
0299: : data[0].length)
0300: + padcol, true);
0301: double[] value = storage.value;
0302: int havecols = cols - padcol;
0303:
0304: for (int i = 0; i < data.length; i++) {
0305: System.arraycopy(data[i], 0, value, i * cols, havecols);
0306: }
0307: }
0308:
0309: /**
0310: * Construct a column vector from a 1D array of doubles. This converts a 1D
0311: * array into a true matrix object. NB. if you want a row vector remember to
0312: * transpose it.
0313: *
0314: * @param data
0315: */
0316: public Matrix(double[] data) {
0317: this (data.length, 1);
0318: System.arraycopy(data, 0, storage.value, 0, data.length);
0319: }
0320:
0321: /**
0322: * Create an identity matrix.
0323: *
0324: * @param n size, both rows and columns
0325: */
0326: static public Matrix identity(int n) {
0327:
0328: int i;
0329:
0330: Matrix ret = new Matrix(n);
0331: for (i = 0; i < n; i++) {
0332: ret.setMval(i, i, 1.0);
0333: }
0334: return ret;
0335: }
0336:
0337: /**
0338: * Create a square matrix of all ones.
0339: *
0340: * @param n the dimension of the required square matrix.
0341: */
0342: static public Matrix ones(int n) {
0343: return (ones(n, n));
0344: }
0345:
0346: /**
0347: * Create a rectangular matrix of all ones.
0348: *
0349: * @param row Number of rows.
0350: * @param col Number of columns.
0351: */
0352: static public Matrix ones(int row, int col) {
0353:
0354: Matrix ret = new Matrix(row, col);
0355: for (int i = 0; i < row; i++) {
0356: for (int j = 0; j < col; j++) {
0357: ret.setMval(i, j, 1.0);
0358: }
0359: }
0360: return ret;
0361: }
0362:
0363: /**
0364: * Create a square matrix of all zeros.
0365: *
0366: * @param n the dimension of the required square matrix.
0367: */
0368: static public Matrix zeros(int n) {
0369: return (zeros(n, n));
0370: }
0371:
0372: /**
0373: * Create a matrix of all zeros.
0374: *
0375: * @param row Number of rows.
0376: * @param col Number of columns.
0377: */
0378: static public Matrix zeros(int row, int col) {
0379:
0380: Matrix ret = new Matrix(row, col);
0381: return ret;
0382: }
0383:
0384: /**
0385: * Clone makes a deep copy of the matrix contents.
0386: */
0387:
0388: public Object clone() {
0389: try {
0390: Matrix res = (Matrix) super .clone();
0391: res.storage = storage;
0392: storage.references++;
0393: // res.value = new double[value.length];
0394: // System.arraycopy(value, 0, res.value, 0, value.length);
0395: return res;
0396: } catch (CloneNotSupportedException e) {
0397: // will never happen
0398: return null;
0399: }
0400: }
0401:
0402: /**
0403: * Perform a deep copy of the matrix.
0404: */
0405: public Matrix copy() {
0406:
0407: Matrix ret = (Matrix) this .clone();
0408: return ret;
0409: }
0410:
0411: /**
0412: * @return
0413: */
0414: public Matrix copyUnshared() {
0415: Matrix ret = (Matrix) this .clone();
0416: ret.aboutToModify();
0417: return ret;
0418: }
0419:
0420: /*
0421: * Basic Matrix Calculations
0422: */
0423:
0424: /**
0425: * Transpose: X'.
0426: *
0427: * @returns A copy of the matrix, transposed.
0428: */
0429: public Matrix transpose() {
0430: int i, j;
0431: // We will create a new matrix which initially shares storage
0432: // the original.
0433: Matrix ret = new Matrix(cols, rows, false);
0434: ret.storage = storage;
0435: storage.references++;
0436: // flip round the indexes.
0437: ret.xstep = ystep;
0438: ret.ystep = xstep;
0439: return ret;
0440: }
0441:
0442: /*
0443: * Transpose the Matrix in place: X' @returns The original matrix, but
0444: * transposed.
0445: */
0446: // AMB NB - for 2x2 matrices and smaller, this is not worth it!
0447: public Matrix updateTranspose() {
0448: int temp = xstep;
0449: xstep = ystep;
0450: ystep = temp;
0451: temp = cols;
0452: cols = rows;
0453: rows = temp;
0454: return this ;
0455: }
0456:
0457: /**
0458: * Linear conversion: a*X + b*Y.
0459: */
0460: public Matrix linearConv(double a, double b, Matrix Y)
0461: throws SizeMismatchException {
0462: // Unoptimised
0463: int i, j;
0464: int Y_rows = Y.rows;
0465: int Y_cols = Y.cols;
0466:
0467: if (rows == Y_rows && cols == Y_cols) {
0468: Matrix ret = new Matrix(rows, cols);
0469:
0470: for (i = 0; i < rows; i++) {
0471: for (j = 0; j < cols; j++) {
0472: ret.setMval(i, j, a * getMval(i, j) + b
0473: * Y.getMval(i, j));
0474: }
0475: }
0476: return ret;
0477: } else {
0478: throw new SizeMismatchException();
0479: }
0480: }
0481:
0482: /**
0483: * Update by Linear conversion: X = a*X + b*Y.
0484: */
0485: public void updateLinearConv(double a, double b, Matrix Y)
0486: throws SizeMismatchException {
0487: // Unoptimised.
0488: int i, j;
0489: int Y_rows = Y.rows;
0490: int Y_cols = Y.cols;
0491:
0492: if (rows != Y_rows || cols != Y_cols) {
0493: throw new SizeMismatchException();
0494: }
0495:
0496: for (i = 0; i < rows; i++) {
0497: for (j = 0; j < cols; j++) {
0498: setMval(i, j, a * getMval(i, j) + b * Y.getMval(i, j));
0499: }
0500: }
0501: }
0502:
0503: /**
0504: * Scalar product: B = aA;
0505: *
0506: * @param a the scalar constant
0507: * @returns a new matrix which is the result of the scalar multiply.
0508: */
0509:
0510: public Matrix multipliedBy(double a) {
0511: Matrix B = copy();
0512: return multipliedByInto(B, a);
0513: }
0514:
0515: /**
0516: * Scalar product: B = aA;
0517: *
0518: * @param a the scalar constant
0519: * @param B the target matrix to be overwritten
0520: * @returns the matrix B which is the result of the scalar multiply.
0521: */
0522:
0523: public Matrix multipliedByInto(Matrix B, double a) {
0524:
0525: int i, j, Aindex = topleft, Bindex = B.topleft;
0526:
0527: B.aboutToModify();
0528:
0529: double[] Aarray = storage.value;
0530: double[] Barray = B.storage.value;
0531:
0532: int Arowspan = (ystep - cols * xstep);
0533: int Browspan = (B.ystep - B.cols * B.xstep);
0534:
0535: for (i = rows; i > 0; i--) {
0536: for (j = cols; j > 0; j--) {
0537: Barray[Bindex] = a * Aarray[Aindex];
0538: Bindex += B.xstep;
0539: Aindex += xstep;
0540: }
0541: Bindex += Browspan;
0542: Aindex += Arowspan;
0543: }
0544: return B;
0545: }
0546:
0547: /**
0548: * Matrix multiplication: C=A*B.
0549: *
0550: * @param B the matrix to multiply by
0551: * @returns a new matrix which is the result of the matrix multiply.
0552: */
0553:
0554: public Matrix multipliedBy(Matrix B) {
0555: Matrix C = new Matrix(rows, B.cols);
0556: return multipliedByInto(C, B);
0557: }
0558:
0559: /**
0560: * Matrix multiplication: C=A*B.
0561: *
0562: * @param C the target matrix to be overwritten by the result
0563: * @param B the matrix to multiply by
0564: * @returns a the matrix C which is the result of the matrix multiply.
0565: */
0566:
0567: public Matrix multipliedByInto(Matrix C, Matrix B) {
0568: int i, j, k; // These can use up the 4 local variable slots in many JVMs.
0569: int B_rows = B.rows, B_cols = B.cols;
0570:
0571: /*
0572: * if (B == C) { throw new UnsupportedOperationException("Cannot multiply
0573: * matrices in place"); }
0574: */
0575:
0576: if (cols != B_rows || C.rows != rows || C.cols != B.cols) {
0577: throw new SizeMismatchException();
0578: }
0579: // this call will produce a matrix that has compacted storage,
0580: // but not necessarily one that is not transposed (i.e. ystep may be
0581: // 1, rather than xstep).
0582: C.aboutToModify();
0583:
0584: double[] Aarray = storage.value;
0585: double[] Barray = B.storage.value;
0586: double[] Carray = C.storage.value;
0587:
0588: // massive strength reduction: instead of a multiply and add
0589: // for each lookup (i.e. 3n^3 integer multiplies and 4n^3 or so adds),
0590: // we set up three registers with the steps required when
0591: // getting to the end of a run on a matrix.
0592: // We then maintain these 3 indices into the 3 matrix arrays, update them
0593: // when required by adding the constant offsets within a row or column.
0594: // When we get to the end of a row or column, we add on the register
0595: // (span) values. Cost is now only 3 int multiplies, and 2n^3 or so
0596: // adds... Of course the n^3 fp multiplies are still all there.
0597:
0598: int Aindex = topleft;
0599: int Bindex = B.topleft;
0600: int Cindex = 0;
0601:
0602: // offset required to move from position off the end of one row
0603: // onto the beginning of the same row of matrix A - this
0604: // is used most frequently in the loop!
0605: int Arowbackfly = -cols * xstep;
0606: // offset required to move from position off the end of one column
0607: // onto the next column of matrix B.
0608: int Bcolspan = (B.xstep - B_rows * B.ystep);
0609: // ditto C - cannot assume it is in good order, since copy()
0610: // may still produce a matrix that is transposed - perhaps think further.
0611: // see aboutToObliterate above, and replace all of these offsets with 1.
0612: int Crowspan = (C.ystep - C.cols * C.xstep);
0613: // i scans down the rows of first matrix
0614: // j scans across the columns of the second matrix
0615: // k is shared between the cols of first matrix and rows of second matrix
0616: // Therefore: for each scan of a row of first matrix, we have a complete
0617: // scan of second matrix.
0618: // Note that the actual values of i,j and k are now irrelevant to the
0619: // calculations - they are only used to time the addition of the constants
0620: // to the indexes. We can therefore make use of the "comparison with 0"
0621: // dodge.
0622: for (i = rows; i > 0; i--) {
0623: for (j = B_cols; j > 0; j--) {
0624: double s = 0;
0625: for (k = cols; k > 0; k--) {
0626: // System.out.println("A: "+Aindex+" B: "+Bindex+" C: "+Cindex);
0627: s += Aarray[Aindex] * Barray[Bindex];
0628: // end for a shared index.
0629: Aindex += xstep; // move the A pointer one to the right.
0630: Bindex += B.ystep; // move the B pointer one down.
0631: }
0632: Carray[Cindex] = s;
0633: // end for a column of the second matrix.
0634: Aindex += Arowbackfly; // Reset the A pointer to the beginning of the
0635: // row.
0636: Bindex += Bcolspan; // Start the second matrix at the next column
0637: Cindex += C.xstep; // Move the output register onto the next element.
0638: }
0639: // end for a row of the first matrix, and hence for all of the second
0640: // matrix.
0641: Bindex = B.topleft; // Start the second matrix at the top left again.
0642: // NB hideous cumulative effects resulting in next line.
0643: Aindex += ystep; // Move the first matrix onto next row, from pos just
0644: // above.
0645: }
0646:
0647: return C;
0648: }
0649:
0650: /**
0651: * Matrix subtraction: C = A-B
0652: */
0653:
0654: public Matrix subtract(Matrix B) {
0655: Matrix C = new Matrix(B.rows, B.cols);
0656: return subtractinto(B, C);
0657: }
0658:
0659: /**
0660: * Matrix subtraction: C = A-B
0661: */
0662:
0663: public Matrix subtractinto(Matrix B, Matrix C) {
0664: int B_rows = B.rows, B_cols = B.cols;
0665:
0666: if (cols != B_cols || rows != B_rows) {
0667: throw new SizeMismatchException();
0668: }
0669: C.aboutToModify();
0670:
0671: int Aindex = topleft;
0672: int Bindex = B.topleft;
0673: int Cindex = C.topleft;
0674:
0675: double[] Aarray = storage.value;
0676: double[] Barray = B.storage.value;
0677: double[] Carray = C.storage.value;
0678:
0679: int Arowspan = (ystep - cols * xstep);
0680: int Browspan = (B.ystep - cols * xstep);
0681: int Crowspan = (C.ystep - cols * xstep);
0682:
0683: for (int i = rows; i > 0; i--) {
0684: for (int j = cols; j > 0; j--) {
0685: Carray[Cindex] = Aarray[Aindex] - Barray[Bindex];
0686: Cindex += C.xstep;
0687: Bindex += B.xstep;
0688: Aindex += xstep;
0689: }
0690: Cindex += Crowspan;
0691: Bindex += Browspan;
0692: Aindex += Arowspan;
0693: }
0694: return C;
0695: }
0696:
0697: public Matrix add(Matrix B) throws SizeMismatchException {
0698: Matrix C = new Matrix(B.rows, B.cols);
0699: return addinto(B, C);
0700: }
0701:
0702: /**
0703: * Stores the addition of this matrix into the first argument into the second
0704: * argument. Either B or A may be equal to C.
0705: * Matrix addition: A+B = C
0706: */
0707:
0708: public Matrix addinto(Matrix B, Matrix C) {
0709: int B_rows = B.rows, B_cols = B.cols;
0710:
0711: if (cols != B_cols || rows != B_rows) {
0712: throw new SizeMismatchException();
0713: }
0714: C.aboutToModify();
0715:
0716: int Aindex = topleft;
0717: int Bindex = B.topleft;
0718: int Cindex = C.topleft;
0719:
0720: double[] Aarray = storage.value;
0721: double[] Barray = B.storage.value;
0722: double[] Carray = C.storage.value;
0723:
0724: int Arowspan = (ystep - cols * xstep);
0725: int Browspan = (B.ystep - cols * xstep);
0726: int Crowspan = (C.ystep - cols * xstep);
0727:
0728: for (int i = rows; i > 0; i--) {
0729: for (int j = cols; j > 0; j--) {
0730: Carray[Cindex] = Aarray[Aindex] + Barray[Bindex];
0731: Cindex += C.xstep;
0732: Bindex += B.xstep;
0733: Aindex += xstep;
0734: }
0735: Cindex += Crowspan;
0736: Bindex += Browspan;
0737: Aindex += Arowspan;
0738: }
0739: return C;
0740: }
0741:
0742: /**
0743: * Matrix left division: X\Y. X\Y is the matrix division of X into Y, which is
0744: * roughly the same as INV(X)*Y , except it is computed in a different way. If
0745: * X is an N-by-N matrix and Y is a column vector with N components, or a
0746: * matrix with several such columns, then B = X\Y is the solution to the
0747: * equation X*B = Y.
0748: */
0749: public Matrix dividedBy(Matrix Y) {
0750: return this .dividedBy(Y, 2 * cols - 1);
0751: }
0752:
0753: /**
0754: * Matrix division with bandwidth: X\Y.
0755: */
0756: public Matrix dividedBy(Matrix Y, int bw) {
0757:
0758: // NB AMB - Optimising this method could prove a real headache!
0759: int i, j, k;
0760:
0761: Matrix L = copy();
0762: Matrix U = Y.copy();
0763:
0764: int m = L.rows;
0765: int n = U.cols;
0766: bw = bw / 2 + 1;
0767:
0768: for (k = 0; k < m; k++) {
0769:
0770: double pp = L.getMval(k, k);
0771:
0772: if (Math.abs(pp) < small) {
0773: throw new UniversalRuntimeException(
0774: "Overflow in Matrix.dividedBy: " + pp);
0775: }
0776:
0777: double w = 1 / pp;
0778:
0779: for (j = k + 1; j < k + bw && j < m; j++) {
0780:
0781: double s = w * L.getMval(k, j);
0782: L.setMval(k, j, s);
0783:
0784: for (i = k + 1; i < k + bw + 1 && i < m; i++) {
0785: s = L.getMval(i, j) - L.getMval(i, k)
0786: * L.getMval(k, j);
0787: L.setMval(i, j, s);
0788: }
0789: }
0790:
0791: for (j = 0; j < n; j++) {
0792:
0793: double s = w * U.getMval(k, j);
0794: U.setMval(k, j, s);
0795:
0796: for (i = k + 1; i < k + bw && i < m; i++) {
0797: s = U.getMval(i, j) - L.getMval(i, k)
0798: * U.getMval(k, j);
0799: U.setMval(i, j, s);
0800: }
0801: }
0802:
0803: }
0804:
0805: Matrix ret = new Matrix(m, n);
0806:
0807: for (i = m - 1; i >= 0; i--) {
0808: for (j = 0; j < n; j++) {
0809:
0810: double w = U.getMval(i, j);
0811:
0812: for (k = i + 1; i < k + bw && k < m; k++) {
0813: w -= L.getMval(i, k) * ret.getMval(k, j);
0814: }
0815:
0816: ret.setMval(i, j, w);
0817: }
0818: }
0819: return ret;
0820: }
0821:
0822: /**
0823: * Cross squared Euclidean distance Matrix: Z_{ij} = ||X_i - Y_j||^2
0824: */
0825: public Matrix crossSqDistance(Matrix Y)
0826: throws SizeMismatchException {
0827: // Unoptimised. In fact, I can't even quite work out what it does!
0828: int i, j, k;
0829: int Y_rows = Y.rows, Y_cols = Y.cols;
0830:
0831: if (cols != Y_cols) {
0832: throw new SizeMismatchException();
0833: }
0834:
0835: Matrix ret = new Matrix(rows, Y_rows);
0836:
0837: for (i = 0; i < rows; i++) {
0838: for (j = 0; j < Y_rows; j++) {
0839: double s = 0;
0840: for (k = 0; k < cols; k++) {
0841: double d = getMval(i, k) - Y.getMval(j, k);
0842: s += d * d;
0843: }
0844: ret.setMval(i, j, s);
0845: }
0846: }
0847: return ret;
0848: }
0849:
0850: public Matrix mapEntries(UnaryFunction func) {
0851:
0852: aboutToModify();
0853:
0854: int i, j, index = topleft;
0855:
0856: double[] array = storage.value;
0857:
0858: int Arowspan = (ystep - cols * xstep);
0859:
0860: for (i = rows; i > 0; i--) {
0861: for (j = cols; j > 0; j--) {
0862: array[index] = func.apply(array[index]);
0863: index += xstep;
0864: }
0865: index += Arowspan;
0866: }
0867: return this ;
0868: }
0869:
0870: /**
0871: * Update by square: (X_{ij})^2
0872: */
0873:
0874: public Matrix updateSqr() {
0875: // Unoptimised
0876: aboutToModify();
0877: for (int i = 0; i < rows; i++) {
0878: for (int j = 0; j < cols; j++) {
0879: double val = getMval(i, j);
0880: setMval(i, j, val * val);
0881: }
0882: }
0883: return this ;
0884: }
0885:
0886: /**
0887: * Create a horizontal-summation vector.
0888: */
0889: public Matrix horizontalSum() {
0890: int i, j, Aindex = topleft;
0891: int vecindex = 0;
0892: int Arowspan = ystep - cols * xstep;
0893: double[] Aarray = storage.value;
0894: Matrix ret = new Matrix(rows, 1);
0895: double[] vecarray = ret.storage.value;
0896:
0897: for (i = rows; i > 0; i--) {
0898: double s = 0;
0899: for (j = cols; j > 0; j--) {
0900: s += Aarray[Aindex];
0901: Aindex += xstep;
0902: }
0903: Aindex += Arowspan;
0904: vecarray[vecindex] = s;
0905: ++vecindex;
0906: }
0907: return ret;
0908: }
0909:
0910: /**
0911: * Create a vertical-summation vector.
0912: */
0913: public Matrix verticalSum() {
0914: int i, j, Aindex = topleft;
0915: int vecindex = 0;
0916: int Acolspan = xstep - rows * ystep;
0917: double[] Aarray = storage.value;
0918: Matrix ret = new Matrix(1, cols);
0919: double[] vecarray = ret.storage.value;
0920:
0921: for (j = cols; j > 0; j--) {
0922: double s = 0;
0923: for (i = rows; i > 0; i--) {
0924: s += Aarray[Aindex];
0925: Aindex += ystep;
0926: }
0927: Aindex += Acolspan;
0928: vecarray[vecindex] = s;
0929: ++vecindex;
0930: }
0931: return ret;
0932: }
0933:
0934: /**
0935: * Summation of all entries.
0936: */
0937: public double sumEntries() {
0938:
0939: int i, j, index = topleft;
0940:
0941: double total = 0;
0942:
0943: double[] array = storage.value;
0944:
0945: int Arowspan = (ystep - cols * xstep);
0946:
0947: for (i = rows; i > 0; i--) {
0948: for (j = cols; j > 0; j--) {
0949: total += array[index];
0950: index += xstep;
0951: }
0952: index += Arowspan;
0953: }
0954: return total;
0955: }
0956:
0957: /**
0958: * Summation of squared entries.
0959: */
0960: public double sumSqrEntries() {
0961:
0962: int i, j, Aindex = topleft;
0963:
0964: double total = 0;
0965:
0966: double[] Aarray = storage.value;
0967:
0968: int Arowspan = (ystep - cols * xstep);
0969:
0970: for (i = rows; i > 0; i--) {
0971: for (j = cols; j > 0; j--) {
0972: double x = Aarray[Aindex];
0973: total += x * x;
0974: Aindex += xstep;
0975: }
0976: Aindex += Arowspan;
0977: }
0978: return total;
0979: }
0980:
0981: /**
0982: * Elementwise multiplication.
0983: */
0984: public Matrix multipliedEntriesBy(Matrix B) {
0985: Matrix C = new Matrix(B.cols, B.rows);
0986: multipliedEntriesByInto(C, B);
0987: return C;
0988: }
0989:
0990: public Matrix multipliedEntriesByInto(Matrix C, Matrix B) {
0991: int i, j;
0992:
0993: if (rows != B.rows || cols != B.cols || cols != C.cols
0994: || rows != C.rows) {
0995: throw new SizeMismatchException();
0996: }
0997: C.aboutToModify();
0998:
0999: int Aindex = topleft;
1000: int Bindex = B.topleft;
1001: int Cindex = C.topleft;
1002:
1003: double[] Aarray = storage.value;
1004: double[] Barray = B.storage.value;
1005: double[] Carray = C.storage.value;
1006:
1007: int Arowspan = (ystep - cols * xstep);
1008: int Browspan = (B.ystep - cols * B.xstep);
1009: int Crowspan = (C.ystep - cols * C.xstep);
1010:
1011: for (i = rows; i > 0; i--) {
1012: for (j = 0; j < cols; j++) {
1013: Carray[Cindex] = Aarray[Aindex] * Barray[Bindex];
1014: Aindex += xstep;
1015: Bindex += B.xstep;
1016: Cindex += C.xstep;
1017: }
1018: Aindex += Arowspan;
1019: Bindex += Browspan;
1020: Cindex += Crowspan;
1021: }
1022: return C;
1023: }
1024:
1025: /**
1026: * X_{ij} = X_{ij}/v_{j}.
1027: */
1028: public Matrix updateScaleRowsByInvVector(Matrix vec) {
1029:
1030: int i, j;
1031:
1032: if (rows != vec.rows || vec.cols > 1) {
1033: throw new SizeMismatchException();
1034: }
1035: aboutToModify();
1036:
1037: double[] Aarray = storage.value;
1038: double[] vecarray = vec.storage.value;
1039: int Aindex = topleft;
1040: int vecindex = vec.topleft;
1041: int Arowspan = (ystep - cols * xstep);
1042:
1043: for (i = rows; i > 0; i--) {
1044: for (j = cols; j > 0; j--) {
1045: Aarray[Aindex] /= vecarray[vecindex];
1046: if (Double.isNaN(Aarray[Aindex]))
1047: Aarray[Aindex] = 0;
1048: Aindex += xstep;
1049: }
1050: vecindex += vec.ystep;
1051: Aindex += Arowspan;
1052: }
1053: return this ;
1054: }
1055:
1056: /**
1057: * X_{ij} = X_{ij}*v_{j}, no summation.
1058: */
1059: // Should really make into and non-into versions of these.
1060: public void updateScaleRowsByColVector(Matrix vec) {
1061:
1062: int i, j;
1063:
1064: if (rows != vec.rows || vec.cols > 1) {
1065: throw new SizeMismatchException();
1066: }
1067: aboutToModify();
1068:
1069: double[] Aarray = storage.value;
1070: double[] vecarray = vec.storage.value;
1071: int Aindex = topleft;
1072: int vecindex = vec.topleft;
1073: int Arowspan = (ystep - cols * xstep);
1074:
1075: for (i = rows; i > 0; i--) {
1076: double multand = vecarray[vecindex];
1077: for (j = cols; j > 0; j--) {
1078: Aarray[Aindex] *= multand;
1079: Aindex += xstep;
1080: }
1081: vecindex += vec.ystep;
1082: Aindex += Arowspan;
1083: }
1084: }
1085:
1086: /**
1087: * Subtract a row vector from each row in this matrix. The number of entries
1088: * in the vector must match the number of columns in the matrix. X_{ij} =
1089: * X_{ij}-v_{j}.
1090: */
1091: public Matrix updateSubtractRowVector(Matrix vec) {
1092: int i, j;
1093: if (vec.cols != cols)
1094: throw new SizeMismatchException("Subtracting " + vec.cols
1095: + " from " + cols);
1096:
1097: aboutToModify();
1098:
1099: int colspan = (xstep - rows * ystep);
1100: int vecindex = vec.topleft;
1101: int arrayindex = topleft;
1102:
1103: double[] array = storage.value;
1104: double[] vecarray = vec.storage.value;
1105:
1106: for (j = cols; j > 0; j--) {
1107: double suband = vecarray[vecindex];
1108:
1109: for (i = rows; i > 0; i--) {
1110: array[arrayindex] -= suband;
1111: arrayindex += ystep;
1112: }
1113: vecindex += vec.xstep;
1114: arrayindex += colspan;
1115: }
1116: return this ;
1117: }
1118:
1119: /**
1120: * If the matrix is a row or column vector, make a diagonal matrix from it.
1121: * Otherwise extract its diagonal as a row vector.
1122: */
1123: // AMB - Hmm... That could be an unfortunate duplication of functions...
1124: public Matrix diagonal() {
1125: // Unoptimised
1126: Matrix ret;
1127: int i, n;
1128:
1129: if (rows == 1) {
1130: ret = new Matrix(cols, cols);
1131: for (i = 0; i < cols; i++) {
1132: ret.setMval(i, i, getMval(0, i));
1133: }
1134: } else if (cols == 1) {
1135: ret = new Matrix(rows, rows);
1136: for (i = 0; i < rows; i++) {
1137: ret.setMval(i, i, getMval(i, 0));
1138: }
1139: } else {
1140: n = Math.min(rows, cols);
1141: ret = new Matrix(1, n);
1142: for (i = 0; i < n; i++) {
1143: ret.setMval(0, i, getMval(i, i));
1144: }
1145: }
1146: return ret;
1147: }
1148:
1149: /**
1150: * Compute Tr(AB) - note this is done in n^2 time, as opposed to the n^3 time
1151: * of multiplying AB and taking trace
1152: */
1153:
1154: public double traceProduct(Matrix B) {
1155:
1156: int i, j, Aindex = topleft, Bindex = B.topleft;
1157:
1158: double total = 0;
1159:
1160: double[] Aarray = storage.value;
1161: double[] Barray = B.storage.value;
1162: int Arowspan = (ystep - cols * xstep);
1163: int Bcolspan = (B.xstep - B.rows * B.ystep);
1164:
1165: for (i = rows; i > 0; i--) {
1166: for (j = cols; j > 0; j--) {
1167: total += Aarray[Aindex] * Barray[Bindex];
1168: Aindex += xstep;
1169: Bindex += ystep;
1170: }
1171: Aindex += Arowspan;
1172: Bindex += Bcolspan;
1173: }
1174: return total;
1175: }
1176:
1177: /**
1178: * Invert the matrix if possible. There is lazy evaluation of the LU
1179: * decomposition. So the first routine to require its presence causes its
1180: * creation.
1181: */
1182: public Matrix inverse() throws SingularException {
1183:
1184: if (LUD == null) {
1185: LUD = new LUDecomp(this );
1186: }
1187:
1188: return (LUD.luinvert());
1189: }
1190:
1191: /**
1192: * Compute the determinant of this matrix if possible. This routine also makes
1193: * use of the LU decomposition. It is party to the lazy evaluation of the LUD
1194: * variable if it has not already been done.
1195: */
1196: public double determinant() throws SingularException {
1197: if (LUD == null) {
1198: LUD = new LUDecomp(this );
1199: }
1200:
1201: return (LUD.ludeterminant());
1202: }
1203:
1204: /**
1205: * Return a copy of the array as a 2D array of doubles. This is a new copy of
1206: * the contents not just a reference to the internals of the matrix object.
1207: */
1208: public double[][] asArray() {
1209: // Unoptimised
1210: double[][] res = new double[rows][cols];
1211:
1212: for (int i = 0; i < rows; i++) {
1213: for (int j = 0; j < cols; j++) {
1214: res[i][j] = getMval(i, j);
1215: }
1216: }
1217: return (res);
1218: }
1219:
1220: /**
1221: * Compute and return the eigen values of this matrix. Since the eigen vectors
1222: * are computed as a side effect the results of both computations are cached
1223: * if the other value is requested.
1224: *
1225: * @returns The matrix of eigen values sorted into increasing magnitude.
1226: */
1227: public Matrix eigenValues() {
1228: if (eigenVal == null)
1229: Jacobi(JacobiIter);
1230:
1231: return (eigenVal);
1232: }
1233:
1234: /**
1235: * Compute and return the eigen vectors of this matrix. Since the eigen values
1236: * are computed as a side effect the results of both computations are cached
1237: * if the other is requested.
1238: *
1239: * @returns The matrix of eigen vectors (as columns of the matrix) with the
1240: * same ordering as that found in the eigen values.
1241: */
1242: public Matrix eigenVectors() {
1243: if (eigenVec == null)
1244: Jacobi(JacobiIter);
1245:
1246: return (eigenVec);
1247: }
1248:
1249: /**
1250: * A direct port of the Numerical Recipes version of the Jacobi algorithim.
1251: * Slightly cleaned up in java it does set the class variables for both the
1252: * Eigen Vectors and Values Matrix variables. It is intended to be applied in
1253: * a single shot mode providing lazy evaluation of either. i.e. when it is
1254: * first called it populates these values for later use.
1255: *
1256: * @param iter is the maximum number of Iterations that should be conducted in
1257: * the primary loop before the system throws a
1258: * JacobiIterationsExhaustedException. This is almost certainly very
1259: * bad if it happens. NR suggests a value of 50 for this parameter.
1260: * @exception NotSquareException the routine was called on an inappropriate
1261: * non-square matrix.
1262: * @exception JacobiIterationsExhaustedException the routine has still failed
1263: * to zero the off diagonals even after the maximum Iteration
1264: * count has been exceded.
1265: */
1266: // NB AMB - that's not lazy evaluation, it's caching isn't it?
1267: // NB2 - this is almost certainly too big to go here... I wonder where
1268: // we can put it...
1269: private void Jacobi(int iter) {
1270:
1271: if (rows != cols) {
1272: throw new NotSquareException();
1273: }
1274:
1275: int nrot = 0, n = rows;
1276:
1277: double sm, thresh, theta, tau, t, s, g, h, c;
1278: double[] b, d, z;
1279: double[][] a, v;
1280:
1281: b = new double[n];
1282: d = new double[n];
1283: z = new double[n];
1284:
1285: v = (Matrix.identity(n)).asArray();
1286: a = this .asArray();
1287:
1288: for (int ip = 0; ip < n; ip++) {
1289:
1290: // Copy the diagonal of this matrix into b and d
1291: // and ensure z is initalised to 0.0.
1292: b[ip] = a[ip][ip];
1293: d[ip] = b[ip];
1294: z[ip] = 0.0;
1295: }
1296:
1297: for (int i = 0; i < iter; i++) {
1298:
1299: // Sum the off-diagonal elements.
1300: sm = 0.0;
1301: for (int ip = 0; ip < n - 1; ip++) {
1302: for (int iq = ip + 1; iq < n; iq++) {
1303: sm += Math.abs(a[ip][iq]);
1304: }
1305: }
1306: Logger.println("Jacobi sweep " + i + " offdiagonal sum "
1307: + sm, Logger.DEBUG_SUBATOMIC);
1308: /*
1309: * for (int q = 0; q < n; ++ q) { Logger.print(a[q][q] + " ",
1310: * Logger.DEBUG_SUBATOMIC); } Logger.println("", Logger.DEBUG_SUBATOMIC);
1311: */
1312: if (sm == 0.0) {
1313: // This was the GOTO 99 and the termiation
1314: // point for the routine. Assumes that the
1315: // off diagonal sum will eventually converge on
1316: // zero.
1317:
1318: // Slight variation from the Numerical recipes
1319: // version here. We sort the Eigen values and
1320: // their associated vectors.
1321:
1322: int k;
1323: double p;
1324:
1325: for (i = 0; i < n - 1; i++) {
1326: k = i;
1327: p = d[i];
1328:
1329: for (int j = i + 1; j < n; j++) {
1330: if (d[j] >= p) {
1331: k = j;
1332: p = d[j];
1333: }
1334: }
1335:
1336: if (k != i) {
1337: d[k] = d[i];
1338: d[i] = p;
1339:
1340: for (int j = 0; j < n; j++) {
1341: p = v[j][i];
1342: v[j][i] = v[j][k];
1343: v[j][k] = p;
1344: }
1345: }
1346: }
1347:
1348: eigenVal = new Matrix(d);
1349: eigenVec = new Matrix(v);
1350: return;
1351:
1352: }
1353:
1354: // On the first 3 Iterations (0,1,2)
1355: if (i < 3)
1356: thresh = 0.2 * sm / (n * n);
1357: else
1358: thresh = 0.0;
1359:
1360: for (int ip = 0; ip < n - 1; ip++) {
1361: for (int iq = ip + 1; iq < n; iq++) {
1362: g = 100.0 * Math.abs(a[ip][iq]);
1363:
1364: // After 4 sweeps, skip the rotation if
1365: // the off-diagonal element is small.
1366:
1367: if (i >= 3
1368: && (Math.abs(d[ip]) + g == Math.abs(d[ip]))
1369: && (Math.abs(d[iq]) + g == Math.abs(d[iq])))
1370: a[ip][iq] = 0.0;
1371: else if (Math.abs(a[ip][iq]) > thresh) {
1372: h = d[iq] - d[ip];
1373: if (Math.abs(h) + g == Math.abs(h))
1374: t = a[ip][iq] / h;
1375: else {
1376: theta = 0.5 * h / a[ip][iq];
1377: t = 1.0 / (Math.abs(theta) + Math
1378: .sqrt(1.0 + (theta * theta)));
1379: if (theta < 0.0)
1380: t = -t;
1381: }
1382: c = 1.0 / Math.sqrt(1 + (t * t));
1383: s = t * c;
1384: tau = s / (1.0 + c);
1385: h = t * a[ip][iq];
1386: z[ip] = z[ip] - h;
1387: z[iq] = z[iq] + h;
1388: d[ip] = d[ip] - h;
1389: d[iq] = d[iq] + h;
1390: a[ip][iq] = 0.0;
1391:
1392: for (int j = 0; j < ip - 1; j++) {
1393: // Case of rotations
1394: // 0 <= j < p-1
1395: g = a[j][ip];
1396: h = a[j][iq];
1397:
1398: a[j][ip] = g - s * (h + g * tau);
1399: a[j][iq] = h + s * (g - h * tau);
1400: }
1401:
1402: for (int j = ip + 1; j < iq - 1; j++) {
1403:
1404: // Case of rotations
1405: // p < j < q
1406:
1407: g = a[ip][j];
1408: h = a[j][iq];
1409:
1410: a[ip][j] = g - s * (h + g * tau);
1411: a[j][iq] = h + s * (g - h * tau);
1412: }
1413:
1414: for (int j = iq + 1; j < n; j++) {
1415:
1416: g = a[ip][j];
1417: h = a[iq][j];
1418:
1419: a[ip][j] = g - s * (h + g * tau);
1420: a[iq][j] = h + s * (g - h * tau);
1421: }
1422:
1423: for (int j = 0; j < n; j++) {
1424:
1425: g = v[j][ip];
1426: h = v[j][iq];
1427:
1428: v[j][ip] = g - s * (h + g * tau);
1429: v[j][iq] = h + s * (g - h * tau);
1430: }
1431:
1432: nrot++;
1433: }
1434: }
1435: }
1436: for (int ip = 0; ip < n; ip++) {
1437: b[ip] = b[ip] + z[ip];
1438: d[ip] = b[ip];
1439: z[ip] = 0.0;
1440: }
1441: }
1442: // We should never get here, the loop should have terminated
1443: // well before we reach this point for any sensible iter value
1444: // Numerical recipes uses a rather arbitrary value of 50 but
1445: // we leave it open.
1446:
1447: throw new JacobiIterationsExhaustedException("After: " + iter
1448: + " Iterations");
1449: }
1450:
1451: /**
1452: * Allow the matrix to be converted to a string representation for debugging
1453: * purposes.
1454: *
1455: * @returns String form of the matrix.
1456: */
1457: public String toString(int limit) {
1458: if (limit > 0 && rows * cols > limit)
1459: return "<" + rows + "x" + cols
1460: + " matrix too big for string representation>";
1461: CharWrap res = new CharWrap();
1462:
1463: for (int i = 0; i < rows; i++) {
1464: if (rows > 10) {
1465: res.appendPad(Integer.toString(i + 1), 3,
1466: CharWrap.PAD_LEFT).append(": ");
1467: }
1468: for (int j = 0; j < cols; j++) {
1469: if (j != 0) {
1470:
1471: res.append(fieldSeparator).appendPad(
1472: decFormat.format(getMval(i, j)), 7,
1473: CharWrap.PAD_RIGHT);
1474: } else {
1475: res.append(decFormat.format(getMval(i, j)));
1476: }
1477: }
1478: if (i < rows) {
1479: res.append("\n");
1480: }
1481: }
1482: return res.toString();
1483: }
1484:
1485: public String toString() {
1486: return toString(250);
1487: }
1488:
1489: public void write(Writer w) throws IOException {
1490: for (int i = 0; i < rows; i++) {
1491: for (int j = 0; j < cols; j++) {
1492: w.write(Double.toString(getMval(i, j)));
1493: w.write(" ");
1494: }
1495: w.write("\n");
1496: }
1497: }
1498:
1499: public void debugsize() {
1500: System.out.println(rows + "x" + cols + " xstep: " + xstep
1501: + " ystep: " + ystep + " topleft: " + topleft + " in "
1502: + storage.value.length + " array");
1503: }
1504:
1505: public void debug() {
1506: System.out.println("xstep: " + xstep + " ystep: " + ystep
1507: + " topleft: " + topleft);
1508: for (int i = 0; i < rows; i++) {
1509: for (int j = 0; j < cols; j++) {
1510: if (j != 0)
1511: System.out.print(fieldSeparator + getMval(i, j));
1512: else
1513: System.out.print(getMval(i, j));
1514: }
1515: System.out.println();
1516: }
1517: }
1518:
1519: /**
1520: * Return the number of rows in the matrix.
1521: *
1522: * @returns n the number of rows in the matrix.
1523: */
1524: public int rows() {
1525: return rows;
1526: }
1527:
1528: /**
1529: * Return the number of columns in the matrix.
1530: *
1531: * @returns n the number of columns in the matrix.
1532: */
1533: public int cols() {
1534: return cols;
1535: }
1536:
1537: /**
1538: * Return a submatrix. The indicies of the matrix run in accordance with
1539: * standard notation from 1,1 in the top left corner to m,n in the bottom
1540: * right corner.
1541: *
1542: * @param y the y coordinate of the top left corner
1543: * @param x the x coordinate of the top left corner
1544: * @param width the number of cells to extract in the x direction
1545: * @param height the number of cells to extract in the y direction.
1546: * @returns the submatrix.
1547: */
1548:
1549: // NB AMB - might be better to swap width and height arguments for consistency
1550: // subject to checking all the code which could reference this method.
1551: public Matrix getSubMatrix(int y, int x, int width, int height) {
1552: --y;
1553: --x; // return the arguments to 0-based versions.
1554: // This method now mysteriously short!
1555: if (x < 0 || y < 0) { // Negative or zero index errors.
1556: throw new ArrayIndexOutOfBoundsException(y + "," + x);
1557: }
1558:
1559: if ((x + width) > cols || (y + height) > rows) {
1560: throw new ArrayIndexOutOfBoundsException(y + "," + x + "+"
1561: + width + "," + height);
1562: }
1563:
1564: Matrix res = new Matrix(height, width, false);
1565: res.storage = storage;
1566: storage.references++;
1567: res.topleft = topleft + x * xstep + y * ystep;
1568: res.xstep = xstep;
1569: res.ystep = ystep;
1570: // Hahaaaargh! Our work here is done!
1571: return res;
1572: }
1573:
1574: /**
1575: * Set a submatrix. The indices of the matrix run in accordance with standard
1576: * notation from 1,1 in the top left corner to m,n in the bottom right corner.
1577: *
1578: * @param y the y coordinate of the top left corner
1579: * @param x the x coordinate of the top left corner
1580: * @param data the matrix to be inserted at this location
1581: */
1582: public void setSubMatrix(int y, int x, Matrix data) {
1583: int i, j;
1584: --y;
1585: --x;
1586: // Unoptimised
1587: if (x < 0 || y < 0 || x >= cols || y >= rows) {
1588: throw new ArrayIndexOutOfBoundsException(y + "," + x);
1589: }
1590:
1591: if ((x + data.cols) > cols || (y + data.rows) > rows) {
1592: throw new ArrayIndexOutOfBoundsException(y + "," + x + "+"
1593: + data.rows + "," + data.cols);
1594: }
1595:
1596: aboutToModify();
1597: double[] Aarray = storage.value;
1598: double[] Barray = data.storage.value;
1599:
1600: int Aindex = x * xstep + y * ystep;
1601: int Bindex = data.topleft;
1602:
1603: int Arowspan = ystep - data.cols * xstep;
1604: int Browspan = data.ystep - data.cols * data.xstep;
1605:
1606: for (i = data.rows; i > 0; i--) {
1607: for (j = data.cols; j > 0; j--) {
1608: Aarray[Aindex] = Barray[Bindex];
1609: Aindex += xstep;
1610: Bindex += data.xstep;
1611: }
1612: Aindex += Arowspan;
1613: Bindex += Browspan;
1614: }
1615: }
1616:
1617: /**
1618: * Get the value in the associated matrix element. The index is in accordance
1619: * with the standard notation 1,1 is top left and m,n is bottom right.
1620: *
1621: * @param y row index
1622: * @param x column index
1623: * @returns the value at this location.
1624: */
1625: public double getValue(int y, int x) {
1626: if (x <= 0 || y <= 0) { // Negative or zero index errors.
1627: throw new ArrayIndexOutOfBoundsException(y + "," + x);
1628: }
1629: if (x > cols || y > rows) { // Outside this matrix range.
1630: throw new ArrayIndexOutOfBoundsException(y + "," + x);
1631: }
1632: return getMval(y - 1, x - 1);
1633: }
1634:
1635: /**
1636: * Set the value in the associated matrix element. The index is in accordance
1637: * with the standard notation 1,1 is top left and m,n is bottom right.
1638: *
1639: * @param y row index
1640: * @param x column index
1641: */
1642: public void setValue(int y, int x, double val) {
1643: if (x <= 0 || y <= 0) { // Negative or zero index errors.
1644: throw new ArrayIndexOutOfBoundsException(y + "," + x);
1645: }
1646: if (x > cols || y > rows) { // Outside this matrix range.
1647: throw new ArrayIndexOutOfBoundsException(y + "," + x);
1648: }
1649: aboutToModify();
1650: setMval(y - 1, x - 1, val);
1651: }
1652:
1653: /**
1654: * Writes the supplied value into both (y,x) and (x,y) entries of the matrix,
1655: * keeping it in a symmetric condition.
1656: *
1657: * @param y index 1
1658: * @param x index 2
1659: * @param val The value to be set
1660: */
1661: public void setSymm(int y, int x, double val) {
1662: aboutToModify();
1663: setMval(y - 1, x - 1, val);
1664: setMval(x - 1, y - 1, val);
1665: }
1666:
1667: public void finalize() {
1668: // System.out.println("Matrix finalized: ");
1669: // debug();
1670: storage.references--;
1671: storage = null; // Cast him off into the blue yonder!
1672: }
1673:
1674: public static void main(String[] argv) {
1675: /*
1676: * double[][] arrA = {{1,0},{0,2}}; Matrix A = new Matrix(arrA); A.debug();
1677: * double[][] arrB = {{2,0},{0,1}}; Matrix B = new Matrix(arrB); Matrix C =
1678: * A.multipliedBy(B); C.debug(); Matrix vec = new Matrix(new double[]
1679: * {1,1}); C.updateSubtractRowVector(vec); C.debug(); Matrix top2 =
1680: * C.getSubMatrix(1,1,1,2); System.out.println("Top 2 entries are ");
1681: * top2.debug(); top2.updateSqr(); top2.debug(); C.debug();
1682: */
1683:
1684: /*
1685: * double[][] arrC = {{0,0},{0,1},{1,1},{1,0}}; Matrix C = new Matrix(arrC);
1686: * double[] arrprob = {0.2, 0.1, 0.2, 0.1}; Matrix prob = new
1687: * Matrix(arrprob); Matrix mean = Probability.Probability.mean(C, prob);
1688: * mean.debug(); Matrix sig = Probability.Probability.covariance(C, prob);
1689: * sig.debug(); Probability.normalPDF n = new Probability.normalPDF(mean,
1690: * sig); Matrix D = n.evaluate(C); D.debug();
1691: */
1692: }
1693:
1694: }
|