/** * Returns the one-norm of matrix <tt>A</tt>, which is the maximum absolute column sum. */ public double norm1(DoubleMatrix2D A) { double max = 0; for (int column = A.columns(); --column >=0; ) { max = Math.max(max, norm1(A.viewColumn(column))); } return max; } /**
/** * Returns the one-norm of matrix <tt>A</tt>, which is the maximum absolute column sum. */ public double norm1(DoubleMatrix2D A) { double max = 0; for (int column = A.columns(); --column >=0; ) { max = Math.max(max, norm1(A.viewColumn(column))); } return max; } /**
/** * Convert discrete data (in yDat) to a matrix of dummy variables (stored in dDat) */ private void makeDummy(){ dDat = factory2D.make(n, lsum); for(int i = 0; i < q; i++){ for(int j = 0; j < l[i]; j++){ DoubleMatrix1D curCol = yDat.viewColumn(i).copy().assign(Functions.equals(j+1)); if(curCol.zSum() == 0) throw new IllegalArgumentException("Discrete data is missing a level: variable " + i + " level " + j); dDat.viewColumn(lcumsum[i]+j).assign(curCol); } } }
public static DoubleMatrix1D flatten(DoubleMatrix2D m){ DoubleMatrix1D[] colArray = new DoubleMatrix1D[m.columns()]; for(int i = 0; i < m.columns(); i++){ colArray[i] = m.viewColumn(i); } return DoubleFactory1D.dense.make(colArray); }
/** * A matrix <tt>A</tt> is <i>diagonally dominant by column</i> if the absolute value of each diagonal element is larger than the sum of the absolute values of the off-diagonal elements in the corresponding column. * <tt>returns true if for all i: abs(A[i,i]) > Sum(abs(A[j,i])); j != i.</tt> * Matrix may but need not be square. * <p> * Note: Ignores tolerance. */ public boolean isDiagonallyDominantByColumn(DoubleMatrix2D A) { cern.jet.math.Functions F = cern.jet.math.Functions.functions; double epsilon = tolerance(); int min = Math.min(A.rows(), A.columns()); for (int i = min; --i >= 0; ) { double diag = Math.abs(A.getQuick(i,i)); diag += diag; if (diag <= A.viewColumn(i).aggregate(F.plus,F.abs)) return false; } return true; } /**
/** * A matrix <tt>A</tt> is <i>diagonally dominant by column</i> if the absolute value of each diagonal element is larger than the sum of the absolute values of the off-diagonal elements in the corresponding column. * <tt>returns true if for all i: abs(A[i,i]) > Sum(abs(A[j,i])); j != i.</tt> * Matrix may but need not be square. * <p> * Note: Ignores tolerance. */ public boolean isDiagonallyDominantByColumn(DoubleMatrix2D A) { cern.jet.math.Functions F = cern.jet.math.Functions.functions; double epsilon = tolerance(); int min = Math.min(A.rows(), A.columns()); for (int i = min; --i >= 0; ) { double diag = Math.abs(A.getQuick(i,i)); diag += diag; if (diag <= A.viewColumn(i).aggregate(F.plus,F.abs)) return false; } return true; } /**
/** * checks if yDat is zero indexed and converts to 1 index. zscores x */ private void fixData(){ double ymin = StatUtils.min(flatten(yDat).toArray()); if(ymin < 0 || ymin > 1) throw new IllegalArgumentException("Discrete data must be either zero or one indexed. Found min index: " + ymin); if(ymin==0){ yDat.assign(Functions.plus(1.0)); } //z-score columns of X for(int i = 0; i < p; i++){ xDat.viewColumn(i).assign(StatUtils.standardizeData(xDat.viewColumn(i).toArray())); } }
/** * Normalizes matrix of p(z|u) such that \forall_z: \sum_u p(z|u) = 1. * * @param pu_z normalized matrix of p(z|u) */ protected void normalizePuz(DoubleMatrix2D pu_z) { for (int z = 0; z < pu_z.columns(); z++) { final DoubleMatrix1D pu_Z = pu_z.viewColumn(z); pu_Z.assign(mult(1 / pu_Z.aggregate(plus, identity))); } }
private void calcWeights(){ weights = factory1D.make(p+q); for(int i = 0; i < p; i++){ weights.set(i, StatUtils.sd(xDat.viewColumn(i).toArray())); } for(int j = 0; j < q; j++){ double curWeight = 0; for(int k = 0; k < l[j] ; k++){ double curp = yDat.viewColumn(j).copy().assign(Functions.equals(k+1)).zSum()/(double) n; curWeight += curp*(1-curp); } weights.set(p+j, Math.sqrt(curWeight)); } }
/** * Normalizes matrix of p(i|z) such that \forall_z: \sum_i p(i|z) = 1. * * @param piz normalized matrix of p(i|z) */ @Override protected void normalizePiz(DoubleMatrix2D piz) { for (int i = 0; i < piz.columns(); i++) { DoubleMatrix1D tmp = piz.viewColumn(i); double norm = tmp.aggregate(Functions.plus, Functions.identity); if (norm != 0.0) { tmp.assign(Functions.mult(1 / norm)); } } } }
/** * Normalizes matrix of p(i|z) such that \forall_z: \sum_i p(i|z) = 1. * * @param piz normalized matrix of p(i|z) */ @Override protected void normalizePiz(DoubleMatrix2D piz) { for (int i = 0; i < piz.columns(); i++) { DoubleMatrix1D tmp = piz.viewColumn(i); double norm = tmp.aggregate(Functions.plus, Functions.identity); if (norm != 0.0) { tmp.assign(Functions.mult(1 / norm)); } } } }
/** Modifies the given matrix square matrix <tt>A</tt> such that it is diagonally dominant by row and column, hence non-singular, hence invertible. For testing purposes only. @param A the square matrix to modify. @throws IllegalArgumentException if <tt>!isSquare(A)</tt>. */ public void generateNonSingular(DoubleMatrix2D A) { checkSquare(A); cern.jet.math.Functions F = cern.jet.math.Functions.functions; int min = Math.min(A.rows(), A.columns()); for (int i = min; --i >= 0; ) { A.setQuick(i,i, 0); } for (int i = min; --i >= 0; ) { double rowSum = A.viewRow(i).aggregate(F.plus,F.abs); double colSum = A.viewColumn(i).aggregate(F.plus,F.abs); A.setQuick(i,i, Math.max(rowSum,colSum) + i+1); } } /**
private static <O> void prepareRR1(int L, DoubleMatrix1D w, DoubleMatrix2D gt, DoubleMatrix2D q, int N, Stream<? extends IdxPref> prefs, DoubleUnaryOperator confidence, double lambda) { int K = w.size(); double[][] x = new double[K + N][K]; double[] y = new double[K + N]; double[] c = new double[K + N]; for (int k = 0; k < K; k++) { gt.viewColumn(k).toArray(x[k]); y[k] = 0.0; c[k] = 1.0; } int[] j = {K}; prefs.forEach(iv -> { q.viewRow(iv.v1).toArray(x[j[0]]); double Cui = confidence.applyAsDouble(iv.v2); y[j[0]] = (Cui * iv.v2) / (Cui - 1); c[j[0]] = Cui - 1; j[0]++; }); doRR1(L, w, x, y, c, lambda); }
/** Modifies the given matrix square matrix <tt>A</tt> such that it is diagonally dominant by row and column, hence non-singular, hence invertible. For testing purposes only. @param A the square matrix to modify. @throws IllegalArgumentException if <tt>!isSquare(A)</tt>. */ public void generateNonSingular(DoubleMatrix2D A) { checkSquare(A); cern.jet.math.Functions F = cern.jet.math.Functions.functions; int min = Math.min(A.rows(), A.columns()); for (int i = min; --i >= 0; ) { A.setQuick(i,i, 0); } for (int i = min; --i >= 0; ) { double rowSum = A.viewRow(i).aggregate(F.plus,F.abs); double colSum = A.viewColumn(i).aggregate(F.plus,F.abs); A.setQuick(i,i, Math.max(rowSum,colSum) + i+1); } } /**
/** Generates and returns the (economy-sized) orthogonal factor <tt>Q</tt>. @return <tt>Q</tt> */ public DoubleMatrix2D getQ () { cern.jet.math.Functions F = cern.jet.math.Functions.functions; DoubleMatrix2D Q = QR.like(); //double[][] Q = X.getArray(); for (int k = n-1; k >= 0; k--) { DoubleMatrix1D QRcolk = QR.viewColumn(k).viewPart(k,m-k); Q.setQuick(k,k, 1); for (int j = k; j < n; j++) { if (QR.getQuick(k,k) != 0) { DoubleMatrix1D Qcolj = Q.viewColumn(j).viewPart(k,m-k); double s = QRcolk.zDotProduct(Qcolj); s = -s / QR.getQuick(k,k); Qcolj.assign(QRcolk, F.plusMult(s)); } } } return Q; } /**
/** Generates and returns the (economy-sized) orthogonal factor <tt>Q</tt>. @return <tt>Q</tt> */ public DoubleMatrix2D getQ () { cern.jet.math.Functions F = cern.jet.math.Functions.functions; DoubleMatrix2D Q = QR.like(); //double[][] Q = X.getArray(); for (int k = n-1; k >= 0; k--) { DoubleMatrix1D QRcolk = QR.viewColumn(k).viewPart(k,m-k); Q.setQuick(k,k, 1); for (int j = k; j < n; j++) { if (QR.getQuick(k,k) != 0) { DoubleMatrix1D Qcolj = Q.viewColumn(j).viewPart(k,m-k); double s = QRcolk.zDotProduct(Qcolj); s = -s / QR.getQuick(k,k); Qcolj.assign(QRcolk, F.plusMult(s)); } } } return Q; } /**
private DoubleMatrix2D scale(DoubleMatrix2D x, boolean scale) { for (int j = 0; j < x.columns(); j++) { DoubleArrayList u = new DoubleArrayList(x.viewColumn(j).toArray()); double mean = Descriptive.mean(u); for (int i = 0; i < x.rows(); i++) { x.set(i, j, x.get(i, j) - mean); } if (scale) { double rms = rms(x.viewColumn(j)); for (int i = 0; i < x.rows(); i++) { x.set(i, j, x.get(i, j) / rms); } } } return x; }
/** * Outer product of two vectors; Sets <tt>A[i,j] = x[i] * y[j]</tt>. * * @param x the first source vector. * @param y the second source vector. * @param A the matrix to hold the results. Set this parameter to <tt>null</tt> to indicate that a new result matrix shall be constructed. * @return A (for convenience only). * @throws IllegalArgumentException if <tt>A.rows() != x.size() || A.columns() != y.size()</tt>. */ public DoubleMatrix2D multOuter(DoubleMatrix1D x, DoubleMatrix1D y, DoubleMatrix2D A) { int rows = x.size(); int columns = y.size(); if (A==null) A = x.like2D(rows,columns); if (A.rows() != rows || A.columns() != columns) throw new IllegalArgumentException(); for (int row = rows; --row >= 0; ) A.viewRow(row).assign(y); for (int column = columns; --column >= 0; ) A.viewColumn(column).assign(x, cern.jet.math.Functions.mult); return A; } /**
/** * Outer product of two vectors; Sets <tt>A[i,j] = x[i] * y[j]</tt>. * * @param x the first source vector. * @param y the second source vector. * @param A the matrix to hold the results. Set this parameter to <tt>null</tt> to indicate that a new result matrix shall be constructed. * @return A (for convenience only). * @throws IllegalArgumentException if <tt>A.rows() != x.size() || A.columns() != y.size()</tt>. */ public DoubleMatrix2D multOuter(DoubleMatrix1D x, DoubleMatrix1D y, DoubleMatrix2D A) { int rows = x.size(); int columns = y.size(); if (A==null) A = x.like2D(rows,columns); if (A.rows() != rows || A.columns() != columns) throw new IllegalArgumentException(); for (int row = rows; --row >= 0; ) A.viewRow(row).assign(y); for (int column = columns; --column >= 0; ) A.viewColumn(column).assign(x, cern.jet.math.Functions.mult); return A; } /**
private static DoubleMatrix2D getGt(final DenseDoubleMatrix2D p, final DenseDoubleMatrix2D q, double lambda) { final int K = p.columns(); DenseDoubleMatrix2D A1 = new DenseDoubleMatrix2D(K, K); q.zMult(q, A1, 1.0, 0.0, true, false); for (int k = 0; k < K; k++) { A1.setQuick(k, k, lambda + A1.getQuick(k, k)); } EigenvalueDecomposition eig = new EigenvalueDecomposition(A1); DoubleMatrix1D d = eig.getRealEigenvalues(); DoubleMatrix2D gt = eig.getV(); for (int k = 0; k < K; k++) { double a = sqrt(d.get(k)); gt.viewColumn(k).assign(x -> a * x); } return gt; }