List of usage examples for org.apache.commons.math.linear RealMatrix setSubMatrix
void setSubMatrix(double[][] subMatrix, int row, int column) throws MatrixIndexException;
row, column using data in the input subMatrix array. From source file:eu.amidst.core.exponentialfamily.EF_Normal_NormalParents2.java
/** * {@inheritDoc}/*from www . j a v a2s . c om*/ */ @Override public void updateNaturalFromMomentParameters() { /* * First step: means and convariances */ CompoundVector globalMomentParam = (CompoundVector) this.momentParameters; double mean_X = globalMomentParam.getXYbaseMatrix().getEntry(0); RealVector mean_Y = globalMomentParam.getTheta_beta0BetaRV(); double cov_XX = globalMomentParam.getcovbaseMatrix().getEntry(0, 0) - mean_X * mean_X; RealMatrix cov_YY = globalMomentParam.getcovbaseMatrix().getSubMatrix(1, nOfParents, 1, nOfParents) .subtract(mean_Y.outerProduct(mean_Y)); RealVector cov_XY = globalMomentParam.getcovbaseMatrix().getSubMatrix(0, 0, 1, nOfParents).getRowVector(0) .subtract(mean_Y.mapMultiply(mean_X)); //RealVector cov_YX = cov_XY; //outerProduct transposes the vector automatically /* * Second step: betas and variance */ RealMatrix cov_YYInverse = new LUDecompositionImpl(cov_YY).getSolver().getInverse(); RealVector beta = cov_YYInverse.preMultiply(cov_XY); double beta_0 = mean_X - beta.dotProduct(mean_Y); double variance = cov_XX - beta.dotProduct(cov_XY); /* * Third step: natural parameters (5 in total) */ /* * 1) theta_0 */ double theta_0 = beta_0 / variance; double[] theta_0array = { theta_0 }; /* * 2) theta_0Theta */ double variance2Inv = 1.0 / (2 * variance); RealVector theta_0Theta = beta.mapMultiply(-beta_0 / variance); ((CompoundVector) this.naturalParameters) .setXYbaseVector(new ArrayRealVector(theta_0array, theta_0Theta.getData())); /* * 3) theta_Minus1 */ double theta_Minus1 = -variance2Inv; /* * 4) theta_beta */ RealVector theta_beta = beta.mapMultiply(variance2Inv); /* * 5) theta_betaBeta */ RealMatrix theta_betaBeta = beta.outerProduct(beta).scalarMultiply(-variance2Inv * 2); /* * Store natural parameters */ RealMatrix natural_XY = new Array2DRowRealMatrix(nOfParents + 1, nOfParents + 1); double[] theta_Minus1array = { theta_Minus1 }; RealVector covXY = new ArrayRealVector(theta_Minus1array, theta_beta.getData()); natural_XY.setColumnVector(0, covXY); natural_XY.setRowVector(0, covXY); natural_XY.setSubMatrix(theta_betaBeta.getData(), 1, 1); ((CompoundVector) this.naturalParameters).setcovbaseVector(natural_XY); }
From source file:name.mjw.cytospade.fcsFile.java
/** * getCompensatedEventList ---/*from w w w . j av a 2 s.com*/ * <p> * Returns the event list compensated by the SPILL matrix. * <p> * * @return array of double arrays containing the events. */ public double[][] getCompensatedEventList() { double[][] events = this.getEventList(); if (events.length != this.getNumChannels()) return events; // Unable to extract the underlying events // Convert the SPILL string to a compensation matrix String compString = this.getSpillString(); if (compString == null) return events; // No compensation, just return the events // Split the compensation string into its values // // The basic structure for SPILL* is: // $SPILLOVER/n,string1,string2,...,f1,f2,f3,f4,.../ String[] compValues = compString.split(","); String[] compNames = null; String[] compData = null; int compDataStart = 0; int n = 0; try { // Try to parse the number of acquisition parameters n = Integer.parseInt(compValues[0]); if (n <= 0 || n > this.parameters) throw new NumberFormatException(); } catch (NumberFormatException nfe) { CyLogger.getLogger().error("Failed to parse parameter count in spill string", nfe); return events; } compNames = Arrays.copyOfRange(compValues, 1, n + 1); // Match names in spill string to columns in parameter lists compDataStart = Arrays.asList(this.channelShortname).indexOf(compNames[0]); if (compDataStart < 0) { CyLogger.getLogger().error("Failed to match channel " + compNames[0] + " to parameter in file"); return events; // Failure match spill string names to channels } for (int i = 0; i < n; i++) { if (!compNames[i].equals(this.channelShortname[compDataStart + i])) { CyLogger.getLogger().error("Spill channel are not continguous parameters in file"); return events; // Spill string columns not in order } } // Extract actual compensation data compData = Arrays.copyOfRange(compValues, n + 1, compValues.length); if (compData.length != (n * n)) return events; /** * Populate the compensation matrix --- The values are stored in * row-major order, i.e., the elements in the first row appear * first. */ double[][] matrix = new double[n][n]; // Loop through the array of compensation values for (int i = 0; i < n; i++) { for (int j = 0; j < n; j++) { try { matrix[i][j] = Double.parseDouble(compData[i * n + j]); } catch (NumberFormatException nfe) { // Set default value If a NumberFormatException occurred matrix[i][j] = 0.0d; } } } // Compute the inverse of the compensation data and then apply // to data matrix (which is column major). Specifically compute // transpose(inverse(<SPILL MATRIX>)) * data RealMatrix comp = (new LUDecompositionImpl(new Array2DRowRealMatrix(matrix))).getSolver().getInverse(); RealMatrix data = new BlockRealMatrix(events); data.setSubMatrix( // Update compensated portion of data matrix comp.transpose() .multiply(data.getSubMatrix(compDataStart, compDataStart + n - 1, 0, this.getEventCount() - 1)) .getData(), compDataStart, 0); return data.getData(); }
From source file:gephi.spade.panel.fcsFile.java
/** * getCompensatedEventList ---//from www.j a v a2 s. c o m * <p> * Returns the event list compensated by the SPILL matrix. * <p> * * @return array of double arrays containing the events. */ public double[][] getCompensatedEventList() { double[][] events = this.getEventList(); if (events.length != this.getNumChannels()) return events; // Unable to extract the underlying events // Convert the SPILL string to a compensation matrix String compString = this.getSpillString(); if (compString == null) return events; // No compensation, just return the events // Split the compensation string into its values // // The basic structure for SPILL* is: // $SPILLOVER/n,string1,string2,...,f1,f2,f3,f4,.../ String[] compValues = compString.split(","); String[] compNames = null; String[] compData = null; int compDataStart = 0; int n = 0; try { // Try to parse the number of acquisition parameters n = Integer.parseInt(compValues[0]); if (n <= 0 || n > this.parameters) throw new NumberFormatException(); } catch (NumberFormatException nfe) { //CyLogger.getLogger().error("Failed to parse parameter count in spill string",nfe); return events; } compNames = Arrays.copyOfRange(compValues, 1, n + 1); // Match names in spill string to columns in parameter lists compDataStart = Arrays.asList(this.channelShortname).indexOf(compNames[0]); if (compDataStart < 0) { //CyLogger.getLogger().error("Failed to match channel "+compNames[0]+" to parameter in file"); return events; // Failure match spill string names to channels } for (int i = 0; i < n; i++) { if (!compNames[i].equals(this.channelShortname[compDataStart + i])) { //CyLogger.getLogger().error("Spill channel are not continguous parameters in file"); return events; // Spill string columns not in order } } // Extract actual compensation data compData = Arrays.copyOfRange(compValues, n + 1, compValues.length); if (compData.length != (n * n)) return events; /** * Populate the compensation matrix --- The values are stored in * row-major order, i.e., the elements in the first row appear * first. */ double[][] matrix = new double[n][n]; // Loop through the array of compensation values for (int i = 0; i < n; i++) { for (int j = 0; j < n; j++) { try { matrix[i][j] = Double.parseDouble(compData[i * n + j]); } catch (NumberFormatException nfe) { // Set default value If a NumberFormatException occurred matrix[i][j] = 0.0d; } } } // Compute the inverse of the compensation data and then apply // to data matrix (which is column major). Specifically compute // transpose(inverse(<SPILL MATRIX>)) * data RealMatrix comp = (new LUDecompositionImpl(new Array2DRowRealMatrix(matrix))).getSolver().getInverse(); RealMatrix data = new BlockRealMatrix(events); data.setSubMatrix( // Update compensated portion of data matrix comp.transpose() .multiply(data.getSubMatrix(compDataStart, compDataStart + n - 1, 0, this.getEventCount() - 1)) .getData(), compDataStart, 0); return data.getData(); }