List of usage examples for org.apache.commons.math3.linear RealMatrix getEntry
double getEntry(int row, int column) throws OutOfRangeException;
From source file:com.google.location.lbs.gnss.gps.pseudorange.UserPositionVelocityWeightedLeastSquare.java
/** * Least square solution to calculate the user position given the navigation message, pseudorange * and accumulated delta range measurements. Also calculates user velocity non-iteratively from * Least square position solution./*from ww w .jav a 2 s . com*/ * * <p>The method fills the user position and velocity in ECEF coordinates and receiver clock * offset in meters and clock offset rate in meters per second. * * <p>One can choose between no smoothing, using the carrier phase measurements (accumulated delta * range) or the doppler measurements (pseudorange rate) for smoothing the pseudorange. The * smoothing is applied only if time has changed below a specific threshold since last invocation. * * <p>Source for least squares: * * <ul> * <li>http://www.u-blox.com/images/downloads/Product_Docs/GPS_Compendium%28GPS-X-02007%29.pdf * page 81 - 85 * <li>Parkinson, B.W., Spilker Jr., J.J.: Global positioning system: theory and applications * page 412 - 414 * </ul> * * <p>Sources for smoothing pseudorange with carrier phase measurements: * * <ul> * <li>Satellite Communications and Navigation Systems book, page 424, * <li>Principles of GNSS, Inertial, and Multisensor Integrated Navigation Systems, page 388, * 389. * </ul> * * <p>The function does not modify the smoothed measurement list {@code * immutableSmoothedSatellitesToReceiverMeasurements} * * @param navMessageProto parameters of the navigation message * @param usefulSatellitesToReceiverMeasurements Map of useful satellite PRN to {@link * GpsMeasurementWithRangeAndUncertainty} containing receiver measurements for computing the * position solution. * @param receiverGPSTowAtReceptionSeconds Receiver estimate of GPS time of week (seconds) * @param receiverGPSWeek Receiver estimate of GPS week (0-1024+) * @param dayOfYear1To366 The day of the year between 1 and 366 * @param positionVelocitySolutionECEF Solution array of the following format: * [0-2] xyz solution of user. * [3] clock bias of user. * [4-6] velocity of user. * [7] clock bias rate of user. * @param positionVelocityUncertaintyEnu Uncertainty of calculated position and velocity solution * in meters and mps local ENU system. Array has the following format: * [0-2] Enu uncertainty of position solution in meters * [3-5] Enu uncertainty of velocity solution in meters per second. * @param pseudorangeResidualMeters The pseudorange residual corrected by subtracting expected * psudorange calculated with the use clock bias of the highest elevation satellites. */ public void calculateUserPositionVelocityLeastSquare(GpsNavMessageProto navMessageProto, List<GpsMeasurementWithRangeAndUncertainty> usefulSatellitesToReceiverMeasurements, double receiverGPSTowAtReceptionSeconds, int receiverGPSWeek, int dayOfYear1To366, double[] positionVelocitySolutionECEF, double[] positionVelocityUncertaintyEnu, double[] pseudorangeResidualMeters) throws Exception { // Use PseudorangeSmoother to smooth the pseudorange according to: Satellite Communications and // Navigation Systems book, page 424 and Principles of GNSS, Inertial, and Multisensor // Integrated Navigation Systems, page 388, 389. double[] deltaPositionMeters; List<GpsMeasurementWithRangeAndUncertainty> immutableSmoothedSatellitesToReceiverMeasurements = pseudorangeSmoother .updatePseudorangeSmoothingResult( Collections.unmodifiableList(usefulSatellitesToReceiverMeasurements)); List<GpsMeasurementWithRangeAndUncertainty> mutableSmoothedSatellitesToReceiverMeasurements = Lists .newArrayList(immutableSmoothedSatellitesToReceiverMeasurements); int numberOfUsefulSatellites = getNumberOfUsefulSatellites(mutableSmoothedSatellitesToReceiverMeasurements); // Least square position solution is supported only if 4 or more satellites visible Preconditions.checkArgument(numberOfUsefulSatellites >= MINIMUM_NUMER_OF_SATELLITES, "At least 4 satellites have to be visible... Only 3D mode is supported..."); boolean repeatLeastSquare = false; SatellitesPositionPseudorangesResidualAndCovarianceMatrix satPosPseudorangeResidualAndWeight; boolean isFirstWLS = true; do { // Calculate satellites' positions, measurement residuals per visible satellite and // weight matrix for the iterative least square boolean doAtmosphericCorrections = false; satPosPseudorangeResidualAndWeight = calculateSatPosAndPseudorangeResidual(navMessageProto, mutableSmoothedSatellitesToReceiverMeasurements, receiverGPSTowAtReceptionSeconds, receiverGPSWeek, dayOfYear1To366, positionVelocitySolutionECEF, doAtmosphericCorrections); // Calculate the geometry matrix according to "Global Positioning System: Theory and // Applications", Parkinson and Spilker page 413 RealMatrix covarianceMatrixM2 = new Array2DRowRealMatrix( satPosPseudorangeResidualAndWeight.covarianceMatrixMetersSquare); geometryMatrix = new Array2DRowRealMatrix(calculateGeometryMatrix( satPosPseudorangeResidualAndWeight.satellitesPositionsMeters, positionVelocitySolutionECEF)); RealMatrix weightedGeometryMatrix; RealMatrix weightMatrixMetersMinus2 = null; // Apply weighted least square only if the covariance matrix is not singular (has a non-zero // determinant), otherwise apply ordinary least square. The reason is to ignore reported // signal to noise ratios by the receiver that can lead to such singularities LUDecomposition ludCovMatrixM2 = new LUDecomposition(covarianceMatrixM2); double det = ludCovMatrixM2.getDeterminant(); if (det <= DOUBLE_ROUND_OFF_TOLERANCE) { // Do not weight the geometry matrix if covariance matrix is singular. weightedGeometryMatrix = geometryMatrix; } else { weightMatrixMetersMinus2 = ludCovMatrixM2.getSolver().getInverse(); RealMatrix hMatrix = calculateHMatrix(weightMatrixMetersMinus2, geometryMatrix); weightedGeometryMatrix = hMatrix.multiply(geometryMatrix.transpose()) .multiply(weightMatrixMetersMinus2); } // Equation 9 page 413 from "Global Positioning System: Theory and Applicaitons", Parkinson // and Spilker deltaPositionMeters = GpsMathOperations.matrixByColVectMultiplication(weightedGeometryMatrix.getData(), satPosPseudorangeResidualAndWeight.pseudorangeResidualsMeters); // Apply corrections to the position estimate positionVelocitySolutionECEF[0] += deltaPositionMeters[0]; positionVelocitySolutionECEF[1] += deltaPositionMeters[1]; positionVelocitySolutionECEF[2] += deltaPositionMeters[2]; positionVelocitySolutionECEF[3] += deltaPositionMeters[3]; // Iterate applying corrections to the position solution until correction is below threshold satPosPseudorangeResidualAndWeight = applyWeightedLeastSquare(navMessageProto, mutableSmoothedSatellitesToReceiverMeasurements, receiverGPSTowAtReceptionSeconds, receiverGPSWeek, dayOfYear1To366, positionVelocitySolutionECEF, deltaPositionMeters, doAtmosphericCorrections, satPosPseudorangeResidualAndWeight, weightMatrixMetersMinus2); // We use the first WLS iteration results and correct them based on the ground truth position // and using a clock error computed from high elevation satellites. The first iteration is // used before satellite with high residuals being removed. if (isFirstWLS && truthLocationForCorrectedResidualComputationEcef != null) { // Snapshot the information needed before high residual satellites are removed System.arraycopy( ResidualCorrectionCalculator.calculateCorrectedResiduals(satPosPseudorangeResidualAndWeight, positionVelocitySolutionECEF.clone(), truthLocationForCorrectedResidualComputationEcef), 0 /*source starting pos*/, pseudorangeResidualMeters, 0 /*destination starting pos*/, GpsNavigationMessageStore.MAX_NUMBER_OF_SATELLITES /*length of elements*/); isFirstWLS = false; } repeatLeastSquare = false; int satsWithResidualBelowThreshold = satPosPseudorangeResidualAndWeight.pseudorangeResidualsMeters.length; // remove satellites that have residuals above RESIDUAL_TO_REPEAT_LEAST_SQUARE_METERS as they // worsen the position solution accuracy. If any satellite is removed, repeat the least square repeatLeastSquare = removeHighResidualSats(mutableSmoothedSatellitesToReceiverMeasurements, repeatLeastSquare, satPosPseudorangeResidualAndWeight, satsWithResidualBelowThreshold); } while (repeatLeastSquare); calculateGeoidMeters = false; // The computed ECEF position will be used next to compute the user velocity. // we calculate and fill in the user velocity solutions based on following equation: // Weight Matrix * GeometryMatrix * User Velocity Vector // = Weight Matrix * deltaPseudoRangeRateWeightedMps // Reference: Pratap Misra and Per Enge // "Global Positioning System: Signals, Measurements, and Performance" Page 218. // Get the number of satellite used in Geometry Matrix numberOfUsefulSatellites = geometryMatrix.getRowDimension(); RealMatrix rangeRateMps = new Array2DRowRealMatrix(numberOfUsefulSatellites, 1); RealMatrix deltaPseudoRangeRateMps = new Array2DRowRealMatrix(numberOfUsefulSatellites, 1); RealMatrix pseudorangeRateWeight = new Array2DRowRealMatrix(numberOfUsefulSatellites, numberOfUsefulSatellites); // Correct the receiver time of week with the estimated receiver clock bias receiverGPSTowAtReceptionSeconds = receiverGPSTowAtReceptionSeconds - positionVelocitySolutionECEF[3] / SPEED_OF_LIGHT_MPS; int measurementCount = 0; // Calculate range rates for (int i = 0; i < GpsNavigationMessageStore.MAX_NUMBER_OF_SATELLITES; i++) { if (mutableSmoothedSatellitesToReceiverMeasurements.get(i) != null) { GpsEphemerisProto ephemeridesProto = getEphemerisForSatellite(navMessageProto, i + 1); double pseudorangeMeasurementMeters = mutableSmoothedSatellitesToReceiverMeasurements .get(i).pseudorangeMeters; GpsTimeOfWeekAndWeekNumber correctedTowAndWeek = calculateCorrectedTransmitTowAndWeek( ephemeridesProto, receiverGPSTowAtReceptionSeconds, receiverGPSWeek, pseudorangeMeasurementMeters); // Calculate satellite velocity PositionAndVelocity satPosECEFMetersVelocityMPS = SatellitePositionCalculator .calculateSatellitePositionAndVelocityFromEphemeris(ephemeridesProto, correctedTowAndWeek.gpsTimeOfWeekSeconds, correctedTowAndWeek.weekNumber, positionVelocitySolutionECEF[0], positionVelocitySolutionECEF[1], positionVelocitySolutionECEF[2]); // Calculate satellite clock error rate double satelliteClockErrorRateMps = SatelliteClockCorrectionCalculator .calculateSatClockCorrErrorRate(ephemeridesProto, correctedTowAndWeek.gpsTimeOfWeekSeconds, correctedTowAndWeek.weekNumber); // Fill in range rates. range rate = satellite velocity (dot product) line-of-sight vector rangeRateMps.setEntry(measurementCount, 0, -1 * (satPosECEFMetersVelocityMPS.velocityXMetersPerSec * geometryMatrix.getEntry(measurementCount, 0) + satPosECEFMetersVelocityMPS.velocityYMetersPerSec * geometryMatrix.getEntry(measurementCount, 1) + satPosECEFMetersVelocityMPS.velocityZMetersPerSec * geometryMatrix.getEntry(measurementCount, 2))); deltaPseudoRangeRateMps.setEntry(measurementCount, 0, mutableSmoothedSatellitesToReceiverMeasurements.get(i).pseudorangeRateMps - rangeRateMps.getEntry(measurementCount, 0) + satelliteClockErrorRateMps - positionVelocitySolutionECEF[7]); // Calculate the velocity weight matrix by using 1 / square(Pseudorangerate Uncertainty) // along the diagonal pseudorangeRateWeight.setEntry(measurementCount, measurementCount, 1 / (mutableSmoothedSatellitesToReceiverMeasurements.get(i).pseudorangeRateUncertaintyMps * mutableSmoothedSatellitesToReceiverMeasurements .get(i).pseudorangeRateUncertaintyMps)); measurementCount++; } } RealMatrix weightedGeoMatrix = pseudorangeRateWeight.multiply(geometryMatrix); RealMatrix deltaPseudoRangeRateWeightedMps = pseudorangeRateWeight.multiply(deltaPseudoRangeRateMps); QRDecomposition qrdWeightedGeoMatrix = new QRDecomposition(weightedGeoMatrix); RealMatrix velocityMps = qrdWeightedGeoMatrix.getSolver().solve(deltaPseudoRangeRateWeightedMps); positionVelocitySolutionECEF[4] = velocityMps.getEntry(0, 0); positionVelocitySolutionECEF[5] = velocityMps.getEntry(1, 0); positionVelocitySolutionECEF[6] = velocityMps.getEntry(2, 0); positionVelocitySolutionECEF[7] = velocityMps.getEntry(3, 0); RealMatrix pseudorangeWeight = new LUDecomposition( new Array2DRowRealMatrix(satPosPseudorangeResidualAndWeight.covarianceMatrixMetersSquare)) .getSolver().getInverse(); // Calculate and store the uncertainties of position and velocity in local ENU system in meters // and meters per second. double[] pvUncertainty = calculatePositionVelocityUncertaintyEnu(pseudorangeRateWeight, pseudorangeWeight, positionVelocitySolutionECEF); System.arraycopy(pvUncertainty, 0 /*source starting pos*/, positionVelocityUncertaintyEnu, 0 /*destination starting pos*/, 6 /*length of elements*/); }
From source file:ffnn.FFNNTubesAI.java
@Override public void buildClassifier(Instances i) throws Exception { Instance temp_instance = null;/*from ww w . j ava 2s . c o m*/ RealMatrix error_output; RealMatrix error_hidden; RealMatrix input_matrix; RealMatrix hidden_matrix; RealMatrix output_matrix; Instances temp_instances; int r = 0; Scanner scan = new Scanner(System.in); output_layer = i.numDistinctValues(i.classIndex()); //3 temp_instances = filterNominalNumeric(i); if (output_layer == 2) { Add filter = new Add(); filter.setAttributeIndex("last"); filter.setAttributeName("dummy"); filter.setInputFormat(temp_instances); temp_instances = Filter.useFilter(temp_instances, filter); // System.out.println(temp_instances); for (int j = 0; j < temp_instances.numInstances(); j++) { if (temp_instances.instance(j).value(temp_instances.numAttributes() - 2) == 0) { temp_instances.instance(j).setValue(temp_instances.numAttributes() - 2, 1); temp_instances.instance(j).setValue(temp_instances.numAttributes() - 1, 0); } else { temp_instances.instance(j).setValue(temp_instances.numAttributes() - 2, 0); temp_instances.instance(j).setValue(temp_instances.numAttributes() - 1, 1); } } } //temp_instances.randomize(temp_instances.getRandomNumberGenerator(1)); //System.out.println(temp_instances); input_layer = temp_instances.numAttributes() - output_layer; //4 hidden_layer = 0; while (hidden_layer < 1) { System.out.print("Hidden layer : "); hidden_layer = scan.nextInt(); } int init_hidden = hidden_layer; error_hidden = new BlockRealMatrix(1, hidden_layer); error_output = new BlockRealMatrix(1, output_layer); input_matrix = new BlockRealMatrix(1, input_layer + 1); //Menambahkan bias buildWeight(input_layer, hidden_layer, output_layer); long last_time = System.nanoTime(); double last_error_rate = 1; double best_error_rate = 1; double last_update = System.nanoTime(); // brp iterasi // for( long itr = 0; last_error_rate > 0.001; ++ itr ){ for (long itr = 0; itr < 50000; ++itr) { if (r == 10) { break; } long time = System.nanoTime(); if (time - last_time > 2000000000) { Evaluation eval = new Evaluation(i); eval.evaluateModel(this, i); double accry = eval.correct() / eval.numInstances(); if (eval.errorRate() < last_error_rate) { last_update = System.nanoTime(); if (eval.errorRate() < best_error_rate) SerializationHelper.write(accry + "-" + time + ".model", this); } if (accry > 0) last_error_rate = eval.errorRate(); // 2 minute without improvement restart if (time - last_update > 30000000000L) { last_update = System.nanoTime(); learning_rate = random() * 0.05; hidden_layer = (int) (10 + floor(random() * 15)); hidden_layer = (int) floor((hidden_layer / 25) * init_hidden); if (hidden_layer == 0) { hidden_layer = 1; } itr = 0; System.out.println("RESTART " + learning_rate + " " + hidden_layer); buildWeight(input_layer, hidden_layer, output_layer); r++; } System.out.println(accry + " " + itr); last_time = time; } for (int j = 0; j < temp_instances.numInstances(); j++) { // foward !! temp_instance = temp_instances.instance(j); for (int k = 0; k < input_layer; k++) { input_matrix.setEntry(0, k, temp_instance.value(k)); } input_matrix.setEntry(0, input_layer, 1.0); // bias hidden_matrix = input_matrix.multiply(weight1); for (int y = 0; y < hidden_layer; ++y) { hidden_matrix.setEntry(0, y, sig(hidden_matrix.getEntry(0, y))); } output_matrix = hidden_matrix.multiply(weight2).add(bias2); for (int y = 0; y < output_layer; ++y) { output_matrix.setEntry(0, y, sig(output_matrix.getEntry(0, y))); } // backward << // error layer 2 double total_err = 0; for (int k = 0; k < output_layer; k++) { double o = output_matrix.getEntry(0, k); double t = temp_instance.value(input_layer + k); double err = o * (1 - o) * (t - o); total_err += err * err; error_output.setEntry(0, k, err); } // back propagation layer 2 for (int y = 0; y < hidden_layer; y++) { for (int x = 0; x < output_layer; ++x) { double wold = weight2.getEntry(y, x); double correction = learning_rate * error_output.getEntry(0, x) * hidden_matrix.getEntry(0, y); weight2.setEntry(y, x, wold + correction); } } for (int x = 0; x < output_layer; ++x) { double correction = learning_rate * error_output.getEntry(0, x); // anggap 1 inputnya bias2.setEntry(0, x, bias2.getEntry(0, x) + correction); } // error layer 1 for (int k = 0; k < hidden_layer; ++k) { double o = hidden_matrix.getEntry(0, k); double t = 0; for (int x = 0; x < output_layer; ++x) { t += error_output.getEntry(0, x) * weight2.getEntry(k, x); } double err = o * (1 - o) * t; error_hidden.setEntry(0, k, err); } // back propagation layer 1 for (int y = 0; y < input_layer + 1; ++y) { for (int x = 0; x < hidden_layer; ++x) { double wold = weight1.getEntry(y, x); double correction = learning_rate * error_hidden.getEntry(0, x) * input_matrix.getEntry(0, y); weight1.setEntry(y, x, wold + correction); } } } } }
From source file:edu.cudenver.bios.power.glmm.NonCentralityDistribution.java
/** * Pre-calculate intermediate matrices, perform setup, etc. *//* ww w . j a v a 2s . com*/ private void initialize(Test test, RealMatrix FEssence, RealMatrix FtFinverse, int perGroupN, FixedRandomMatrix CFixedRand, RealMatrix U, RealMatrix thetaNull, RealMatrix beta, RealMatrix sigmaError, RealMatrix sigmaG, boolean exact) throws PowerException { debug("entering initialize"); // reset member variables this.T1 = null; this.FT1 = null; this.S = null; this.mzSq = null; this.H0 = 0; this.sStar = 0; // cache inputs this.test = test; this.FEssence = FEssence; this.FtFinverse = FtFinverse; this.perGroupN = perGroupN; this.CFixedRand = CFixedRand; this.U = U; this.thetaNull = thetaNull; this.beta = beta; this.sigmaError = sigmaError; this.sigmaG = sigmaG; // calculate intermediate matrices // RealMatrix FEssence = params.getDesignEssence().getFullDesignMatrixFixed(); // TODO: do we ever get here with values that can cause integer overflow, // and if so, does it matter? this.N = (double) FEssence.getRowDimension() * perGroupN; this.exact = exact; try { // TODO: need to calculate H0, need to adjust H1 for Unirep // get design matrix for fixed parameters only qF = FEssence.getColumnDimension(); // a = CFixedRand.getCombinedMatrix().getRowDimension(); // get fixed contrasts RealMatrix Cfixed = CFixedRand.getFixedMatrix(); RealMatrix CGaussian = CFixedRand.getRandomMatrix(); // build intermediate terms h1, S if (FtFinverse == null) { FtFinverse = new LUDecomposition(FEssence.transpose().multiply(FEssence)).getSolver().getInverse(); debug("FEssence", FEssence); debug("FtFinverse = (FEssence transpose * FEssence) inverse", FtFinverse); } else { debug("FtFinverse", FtFinverse); } RealMatrix PPt = Cfixed.multiply(FtFinverse.scalarMultiply(1 / (double) perGroupN)) .multiply(Cfixed.transpose()); debug("Cfixed", Cfixed); debug("n = " + perGroupN); debug("PPt = Cfixed * FtF inverse * (1/n) * Cfixed transpose", PPt); T1 = forceSymmetric(new LUDecomposition(PPt).getSolver().getInverse()); debug("T1 = PPt inverse", T1); FT1 = new CholeskyDecomposition(T1).getL(); debug("FT1 = Cholesky decomposition (L) of T1", FT1); // calculate theta difference // RealMatrix thetaNull = params.getTheta(); RealMatrix C = CFixedRand.getCombinedMatrix(); // RealMatrix beta = params.getScaledBeta(); // RealMatrix U = params.getWithinSubjectContrast(); // thetaHat = C * beta * U RealMatrix thetaHat = C.multiply(beta.multiply(U)); debug("C", C); debug("beta", beta); debug("U", U); debug("thetaHat = C * beta * U", thetaHat); // thetaDiff = thetaHat - thetaNull RealMatrix thetaDiff = thetaHat.subtract(thetaNull); debug("thetaNull", thetaNull); debug("thetaDiff = thetaHat - thetaNull", thetaDiff); // TODO: specific to HLT or UNIREP RealMatrix sigmaStarInverse = getSigmaStarInverse(U, sigmaError, test); debug("sigmaStarInverse", sigmaStarInverse); RealMatrix H1matrix = thetaDiff.transpose().multiply(T1).multiply(thetaDiff).multiply(sigmaStarInverse); debug("H1matrix = thetaDiff transpose * T1 * thetaDiff * sigmaStarInverse", H1matrix); H1 = H1matrix.getTrace(); debug("H1 = " + H1); // Matrix which represents the non-centrality parameter as a linear combination of chi-squared r.v.'s. S = FT1.transpose().multiply(thetaDiff).multiply(sigmaStarInverse).multiply(thetaDiff.transpose()) .multiply(FT1).scalarMultiply(1 / H1); debug("S = FT1 transpose * thetaDiff * sigmaStar inverse * thetaDiff transpose * FT1 * (1/H1)", S); // We use the S matrix to generate the F-critical, numerical df's, and denominator df's // for a central F distribution. The resulting F distribution is used as an approximation // for the distribution of the non-centrality parameter. // See formulas 18-21 and A8,A10 from Glueck & Muller (2003) for details. EigenDecomposition sEigenDecomp = new EigenDecomposition(S); sEigenValues = sEigenDecomp.getRealEigenvalues(); // calculate H0 if (sEigenValues.length > 0) H0 = H1 * (1 - sEigenValues[0]); if (H0 <= 0) H0 = 0; // count the # of positive eigen values for (double value : sEigenValues) { if (value > 0) sStar++; } // TODO: throw error if sStar is <= 0 // TODO: NO: throw error if sStar != sEigenValues.length instead??? double stddevG = Math.sqrt(sigmaG.getEntry(0, 0)); RealMatrix svec = sEigenDecomp.getVT(); mzSq = svec.multiply(FT1.transpose()).multiply(CGaussian).scalarMultiply(1 / stddevG); for (int i = 0; i < mzSq.getRowDimension(); i++) { for (int j = 0; j < mzSq.getColumnDimension(); j++) { double entry = mzSq.getEntry(i, j); mzSq.setEntry(i, j, entry * entry); // TODO: is there an apache function to do this? } } debug("exiting initialize normally"); } catch (RuntimeException e) { LOGGER.warn("exiting initialize abnormally", e); throw new PowerException(e.getMessage(), PowerErrorEnum.INVALID_DISTRIBUTION_NONCENTRALITY_PARAMETER); } }
From source file:org.akvo.caddisfly.sensor.colorimetry.strip.calibration.CalibrationCard.java
@NonNull private static Mat doIlluminationCorrection(@NonNull Mat imgLab, @NonNull CalibrationData calData) { // create HLS image for homogeneous illumination calibration int pHeight = imgLab.rows(); int pWidth = imgLab.cols(); RealMatrix points = createWhitePointMatrix(imgLab, calData); // create coefficient matrix for all three variables L,A,B // the model for all three is y = ax + bx^2 + cy + dy^2 + exy + f // 6th row is the constant 1 RealMatrix coefficient = new Array2DRowRealMatrix(points.getRowDimension(), 6); coefficient.setColumnMatrix(0, points.getColumnMatrix(0)); coefficient.setColumnMatrix(2, points.getColumnMatrix(1)); //create constant, x^2, y^2 and xy terms for (int i = 0; i < points.getRowDimension(); i++) { coefficient.setEntry(i, 1, Math.pow(coefficient.getEntry(i, 0), 2)); // x^2 coefficient.setEntry(i, 3, Math.pow(coefficient.getEntry(i, 2), 2)); // y^2 coefficient.setEntry(i, 4, coefficient.getEntry(i, 0) * coefficient.getEntry(i, 2)); // xy coefficient.setEntry(i, 5, 1d); // constant = 1 }//from w w w .j a v a 2s . c om // create vectors RealVector L = points.getColumnVector(2); RealVector A = points.getColumnVector(3); RealVector B = points.getColumnVector(4); // solve the least squares problem for all three variables DecompositionSolver solver = new SingularValueDecomposition(coefficient).getSolver(); RealVector solutionL = solver.solve(L); RealVector solutionA = solver.solve(A); RealVector solutionB = solver.solve(B); // get individual coefficients float La = (float) solutionL.getEntry(0); float Lb = (float) solutionL.getEntry(1); float Lc = (float) solutionL.getEntry(2); float Ld = (float) solutionL.getEntry(3); float Le = (float) solutionL.getEntry(4); float Lf = (float) solutionL.getEntry(5); float Aa = (float) solutionA.getEntry(0); float Ab = (float) solutionA.getEntry(1); float Ac = (float) solutionA.getEntry(2); float Ad = (float) solutionA.getEntry(3); float Ae = (float) solutionA.getEntry(4); float Af = (float) solutionA.getEntry(5); float Ba = (float) solutionB.getEntry(0); float Bb = (float) solutionB.getEntry(1); float Bc = (float) solutionB.getEntry(2); float Bd = (float) solutionB.getEntry(3); float Be = (float) solutionB.getEntry(4); float Bf = (float) solutionB.getEntry(5); // compute mean (the luminosity value of the plane in the middle of the image) float L_mean = (float) (0.5 * La * pWidth + 0.5 * Lc * pHeight + Lb * pWidth * pWidth / 3.0 + Ld * pHeight * pHeight / 3.0 + Le * 0.25 * pHeight * pWidth + Lf); float A_mean = (float) (0.5 * Aa * pWidth + 0.5 * Ac * pHeight + Ab * pWidth * pWidth / 3.0 + Ad * pHeight * pHeight / 3.0 + Ae * 0.25 * pHeight * pWidth + Af); float B_mean = (float) (0.5 * Ba * pWidth + 0.5 * Bc * pHeight + Bb * pWidth * pWidth / 3.0 + Bd * pHeight * pHeight / 3.0 + Be * 0.25 * pHeight * pWidth + Bf); // Correct image // we do this per row. We tried to do it in one block, but there is no speed difference. byte[] temp = new byte[imgLab.cols() * imgLab.channels()]; int valL, valA, valB; int ii, ii3; float iiSq, iSq; int imgCols = imgLab.cols(); int imgRows = imgLab.rows(); // use lookup tables to speed up computation // create lookup tables float[] L_aii = new float[imgCols]; float[] L_biiSq = new float[imgCols]; float[] A_aii = new float[imgCols]; float[] A_biiSq = new float[imgCols]; float[] B_aii = new float[imgCols]; float[] B_biiSq = new float[imgCols]; float[] Lci = new float[imgRows]; float[] LdiSq = new float[imgRows]; float[] Aci = new float[imgRows]; float[] AdiSq = new float[imgRows]; float[] Bci = new float[imgRows]; float[] BdiSq = new float[imgRows]; for (ii = 0; ii < imgCols; ii++) { iiSq = ii * ii; L_aii[ii] = La * ii; L_biiSq[ii] = Lb * iiSq; A_aii[ii] = Aa * ii; A_biiSq[ii] = Ab * iiSq; B_aii[ii] = Ba * ii; B_biiSq[ii] = Bb * iiSq; } for (int i = 0; i < imgRows; i++) { iSq = i * i; Lci[i] = Lc * i; LdiSq[i] = Ld * iSq; Aci[i] = Ac * i; AdiSq[i] = Ad * iSq; Bci[i] = Bc * i; BdiSq[i] = Bd * iSq; } // We can also improve the performance of the i,ii term, if we want, but it won't make much difference. for (int i = 0; i < imgRows; i++) { // y imgLab.get(i, 0, temp); ii3 = 0; for (ii = 0; ii < imgCols; ii++) { //x valL = capValue( Math.round((temp[ii3] & 0xFF) - (L_aii[ii] + L_biiSq[ii] + Lci[i] + LdiSq[i] + Le * i * ii + Lf) + L_mean), 0, 255); valA = capValue( Math.round((temp[ii3 + 1] & 0xFF) - (A_aii[ii] + A_biiSq[ii] + Aci[i] + AdiSq[i] + Ae * i * ii + Af) + A_mean), 0, 255); valB = capValue( Math.round((temp[ii3 + 2] & 0xFF) - (B_aii[ii] + B_biiSq[ii] + Bci[i] + BdiSq[i] + Be * i * ii + Bf) + B_mean), 0, 255); temp[ii3] = (byte) valL; temp[ii3 + 1] = (byte) valA; temp[ii3 + 2] = (byte) valB; ii3 += 3; } imgLab.put(i, 0, temp); } return imgLab; }
From source file:org.akvo.caddisfly.sensor.colorimetry.strip.calibration.CalibrationCard.java
@NonNull private static Mat do1D_3DCorrection(@NonNull Mat imgMat, @Nullable CalibrationData calData) throws CalibrationException { if (calData == null) { throw new CalibrationException("no calibration data."); }//from ww w.j a va 2s . co m final WeightedObservedPoints obsL = new WeightedObservedPoints(); final WeightedObservedPoints obsA = new WeightedObservedPoints(); final WeightedObservedPoints obsB = new WeightedObservedPoints(); Map<String, double[]> calResultIllumination = new HashMap<>(); // iterate over all patches try { for (String label : calData.getCalValues().keySet()) { CalibrationData.CalValue cal = calData.getCalValues().get(label); CalibrationData.Location loc = calData.getLocations().get(label); float[] LAB_color = measurePatch(imgMat, loc.x, loc.y, calData); // measure patch color obsL.add(LAB_color[0], cal.getL()); obsA.add(LAB_color[1], cal.getA()); obsB.add(LAB_color[2], cal.getB()); calResultIllumination.put(label, new double[] { LAB_color[0], LAB_color[1], LAB_color[2] }); } } catch (Exception e) { throw new CalibrationException("1D calibration: error iterating over all patches.", e); } // Instantiate a second-degree polynomial fitter. final PolynomialCurveFitter fitter = PolynomialCurveFitter.create(2); // Retrieve fitted parameters (coefficients of the polynomial function). // order of coefficients is (c + bx + ax^2), so [c,b,a] try { final double[] coefficientL = fitter.fit(obsL.toList()); final double[] coefficientA = fitter.fit(obsA.toList()); final double[] coefficientB = fitter.fit(obsB.toList()); double[] valIllumination; double L_orig, A_orig, B_orig, L_new, A_new, B_new; // transform patch values using the 1d calibration results Map<String, double[]> calResult1D = new HashMap<>(); for (String label : calData.getCalValues().keySet()) { valIllumination = calResultIllumination.get(label); L_orig = valIllumination[0]; A_orig = valIllumination[1]; B_orig = valIllumination[2]; L_new = coefficientL[2] * L_orig * L_orig + coefficientL[1] * L_orig + coefficientL[0]; A_new = coefficientA[2] * A_orig * A_orig + coefficientA[1] * A_orig + coefficientA[0]; B_new = coefficientB[2] * B_orig * B_orig + coefficientB[1] * B_orig + coefficientB[0]; calResult1D.put(label, new double[] { L_new, A_new, B_new }); } // use the 1D calibration result for the second calibration step // Following http://docs.scipy.org/doc/scipy/reference/tutorial/linalg.html#solving-linear-least-squares-problems-and-pseudo-inverses // we will solve P = M x int total = calData.getLocations().keySet().size(); RealMatrix coefficient = new Array2DRowRealMatrix(total, 3); RealMatrix cal = new Array2DRowRealMatrix(total, 3); int index = 0; // create coefficient and calibration vectors for (String label : calData.getCalValues().keySet()) { CalibrationData.CalValue calv = calData.getCalValues().get(label); double[] cal1dResult = calResult1D.get(label); coefficient.setEntry(index, 0, cal1dResult[0]); coefficient.setEntry(index, 1, cal1dResult[1]); coefficient.setEntry(index, 2, cal1dResult[2]); cal.setEntry(index, 0, calv.getL()); cal.setEntry(index, 1, calv.getA()); cal.setEntry(index, 2, calv.getB()); index++; } DecompositionSolver solver = new SingularValueDecomposition(coefficient).getSolver(); RealMatrix sol = solver.solve(cal); float a_L, b_L, c_L, a_A, b_A, c_A, a_B, b_B, c_B; a_L = (float) sol.getEntry(0, 0); b_L = (float) sol.getEntry(1, 0); c_L = (float) sol.getEntry(2, 0); a_A = (float) sol.getEntry(0, 1); b_A = (float) sol.getEntry(1, 1); c_A = (float) sol.getEntry(2, 1); a_B = (float) sol.getEntry(0, 2); b_B = (float) sol.getEntry(1, 2); c_B = (float) sol.getEntry(2, 2); //use the solution to correct the image double L_temp, A_temp, B_temp, L_mid, A_mid, B_mid; int L_fin, A_fin, B_fin; int ii3; byte[] temp = new byte[imgMat.cols() * imgMat.channels()]; for (int i = 0; i < imgMat.rows(); i++) { // y imgMat.get(i, 0, temp); ii3 = 0; for (int ii = 0; ii < imgMat.cols(); ii++) { //x L_temp = temp[ii3] & 0xFF; A_temp = temp[ii3 + 1] & 0xFF; B_temp = temp[ii3 + 2] & 0xFF; L_mid = coefficientL[2] * L_temp * L_temp + coefficientL[1] * L_temp + coefficientL[0]; A_mid = coefficientA[2] * A_temp * A_temp + coefficientA[1] * A_temp + coefficientA[0]; B_mid = coefficientB[2] * B_temp * B_temp + coefficientB[1] * B_temp + coefficientB[0]; L_fin = (int) Math.round(a_L * L_mid + b_L * A_mid + c_L * B_mid); A_fin = (int) Math.round(a_A * L_mid + b_A * A_mid + c_A * B_mid); B_fin = (int) Math.round(a_B * L_mid + b_B * A_mid + c_B * B_mid); // cap values L_fin = capValue(L_fin, 0, 255); A_fin = capValue(A_fin, 0, 255); B_fin = capValue(B_fin, 0, 255); temp[ii3] = (byte) L_fin; temp[ii3 + 1] = (byte) A_fin; temp[ii3 + 2] = (byte) B_fin; ii3 += 3; } imgMat.put(i, 0, temp); } return imgMat; } catch (Exception e) { throw new CalibrationException("error while performing calibration: ", e); } }
From source file:org.cirdles.calamari.algorithms.WeightedMeanCalculators.java
public static WeightedLinearCorrResults weightedLinearCorr(double[] y, double[] x, double[][] sigmaRhoY) { WeightedLinearCorrResults weightedLinearCorrResults = new WeightedLinearCorrResults(); RealMatrix omega = new BlockRealMatrix(convertCorrelationsToCovariances(sigmaRhoY)); RealMatrix invOmega = MatrixUtils.inverse(omega); int n = y.length; double mX = 0; double pX = 0; double pY = 0; double pXY = 0; double w = 0; for (int i = 0; i < n; i++) { for (int j = 0; j < n; j++) { double invOm = invOmega.getEntry(i, j); w += invOm;//w w w .ja v a 2s .c o m pX += (invOm * (x[i] + x[j])); pY += (invOm * (y[i] + y[j])); pXY += (invOm * (((x[i] * y[j]) + (x[j] * y[i])))); mX += (invOm * x[i] * x[j]); } } double slope = ((2 * pXY * w) - (pX * pY)) / ((4 * mX * w) - (pX * pX)); double intercept = (pY - (slope * pX)) / (2 * w); RealMatrix fischer = new BlockRealMatrix(new double[][] { { mX, pX / 2.0 }, { pX / 2.0, w } }); RealMatrix fischerInv = MatrixUtils.inverse(fischer); double slopeSig = StrictMath.sqrt(fischerInv.getEntry(0, 0)); double interceptSig = StrictMath.sqrt(fischerInv.getEntry(1, 1)); double slopeInterceptCov = fischerInv.getEntry(0, 1); RealMatrix resid = new BlockRealMatrix(n, 1); for (int i = 0; i < n; i++) { resid.setEntry(i, 0, y[i] - (slope * x[i]) - intercept); } RealMatrix residT = resid.transpose(); RealMatrix mM = residT.multiply(invOmega).multiply(resid); double sumSqWtdResids = mM.getEntry(0, 0); double mswd = sumSqWtdResids / (n - 2); // http://commons.apache.org/proper/commons-math/apidocs/org/apache/commons/math3/distribution/FDistribution.html FDistribution fdist = new org.apache.commons.math3.distribution.FDistribution((n - 2), 1E9); double prob = 1.0 - fdist.cumulativeProbability(mswd); weightedLinearCorrResults.setBad(false); weightedLinearCorrResults.setSlope(slope); weightedLinearCorrResults.setIntercept(intercept); weightedLinearCorrResults.setSlopeSig(slopeSig); weightedLinearCorrResults.setInterceptSig(interceptSig); weightedLinearCorrResults.setSlopeInterceptCov(slopeInterceptCov); weightedLinearCorrResults.setMswd(mswd); weightedLinearCorrResults.setProb(prob); return weightedLinearCorrResults; }
From source file:org.cirdles.calamari.algorithms.WeightedMeanCalculators.java
public static WtdAvCorrResults wtdAvCorr(double[] values, double[][] varCov) { // assume varCov is variance-covariance matrix (i.e. SigRho = false) WtdAvCorrResults results = new WtdAvCorrResults(); int n = varCov.length; RealMatrix omegaInv = new BlockRealMatrix(varCov); RealMatrix omega = MatrixUtils.inverse(omegaInv); double numer = 0.0; double denom = 0.0; for (int i = 0; i < n; i++) { for (int j = 0; j < n; j++) { numer += (values[i] + values[j]) * omega.getEntry(i, j); denom += omega.getEntry(i, j); }/*from ww w . j a v a2 s . c o m*/ } // test denom if (denom > 0.0) { double meanVal = numer / denom / 2.0; double meanValSigma = StrictMath.sqrt(1.0 / denom); double[][] unwtdResidsArray = new double[n][1]; for (int i = 0; i < n; i++) { unwtdResidsArray[i][0] = values[i] - meanVal; } RealMatrix unwtdResids = new BlockRealMatrix(unwtdResidsArray); RealMatrix transUnwtdResids = unwtdResids.transpose(); RealMatrix product = transUnwtdResids.multiply(omega); RealMatrix sumWtdResids = product.multiply(unwtdResids); double mswd = 0.0; double prob = 0.0; if (n > 1) { mswd = sumWtdResids.getEntry(0, 0) / (n - 1); // http://commons.apache.org/proper/commons-math/apidocs/org/apache/commons/math3/distribution/FDistribution.html FDistribution fdist = new org.apache.commons.math3.distribution.FDistribution((n - 1), 1E9); prob = 1.0 - fdist.cumulativeProbability(mswd); } results.setBad(false); results.setMeanVal(meanVal); results.setSigmaMeanVal(meanValSigma); results.setMswd(mswd); results.setProb(prob); } return results; }
From source file:org.cirdles.squid.algorithms.weightedMeans.WeightedMeanCalculators.java
/** * * @param y//from w w w . jav a 2 s .co m * @param x * @param sigmaRhoY * @return */ public static WeightedLinearCorrResults weightedLinearCorr(double[] y, double[] x, double[][] sigmaRhoY) { WeightedLinearCorrResults weightedLinearCorrResults = new WeightedLinearCorrResults(); RealMatrix omega = new BlockRealMatrix(convertCorrelationsToCovariances(sigmaRhoY)); RealMatrix invOmega = MatrixUtils.inverse(omega); int n = y.length; double mX = 0; double pX = 0; double pY = 0; double pXY = 0; double w = 0; for (int i = 0; i < n; i++) { for (int j = 0; j < n; j++) { double invOm = invOmega.getEntry(i, j); w += invOm; pX += (invOm * (x[i] + x[j])); pY += (invOm * (y[i] + y[j])); pXY += (invOm * (((x[i] * y[j]) + (x[j] * y[i])))); mX += (invOm * x[i] * x[j]); } } double slope = ((2 * pXY * w) - (pX * pY)) / ((4 * mX * w) - (pX * pX)); double intercept = (pY - (slope * pX)) / (2 * w); RealMatrix fischer = new BlockRealMatrix(new double[][] { { mX, pX / 2.0 }, { pX / 2.0, w } }); RealMatrix fischerInv = MatrixUtils.inverse(fischer); double slopeSig = Math.sqrt(fischerInv.getEntry(0, 0)); double interceptSig = Math.sqrt(fischerInv.getEntry(1, 1)); double slopeInterceptCov = fischerInv.getEntry(0, 1); RealMatrix resid = new BlockRealMatrix(n, 1); for (int i = 0; i < n; i++) { resid.setEntry(i, 0, y[i] - (slope * x[i]) - intercept); } RealMatrix residT = resid.transpose(); RealMatrix mM = residT.multiply(invOmega).multiply(resid); double sumSqWtdResids = mM.getEntry(0, 0); double mswd = sumSqWtdResids / (n - 2); // http://commons.apache.org/proper/commons-math/apidocs/org/apache/commons/math3/distribution/FDistribution.html FDistribution fdist = new org.apache.commons.math3.distribution.FDistribution((n - 2), 1E9); double prob = 1.0 - fdist.cumulativeProbability(mswd); weightedLinearCorrResults.setBad(false); weightedLinearCorrResults.setSlope(slope); weightedLinearCorrResults.setIntercept(intercept); weightedLinearCorrResults.setSlopeSig(slopeSig); weightedLinearCorrResults.setInterceptSig(interceptSig); weightedLinearCorrResults.setSlopeInterceptCov(slopeInterceptCov); weightedLinearCorrResults.setMswd(mswd); weightedLinearCorrResults.setProb(prob); return weightedLinearCorrResults; }
From source file:org.cirdles.squid.algorithms.weightedMeans.WeightedMeanCalculators.java
/** * * @param values/*from ww w. j av a 2 s.c om*/ * @param varCov * @return */ public static WtdAvCorrResults wtdAvCorr(double[] values, double[][] varCov) { // assume varCov is variance-covariance matrix (i.e. SigRho = false) WtdAvCorrResults results = new WtdAvCorrResults(); int n = varCov.length; RealMatrix omegaInv = new BlockRealMatrix(varCov); RealMatrix omega = MatrixUtils.inverse(omegaInv); double numer = 0.0; double denom = 0.0; for (int i = 0; i < n; i++) { for (int j = 0; j < n; j++) { numer += (values[i] + values[j]) * omega.getEntry(i, j); denom += omega.getEntry(i, j); } } // test denom if (denom > 0.0) { double meanVal = numer / denom / 2.0; double meanValSigma = Math.sqrt(1.0 / denom); double[][] unwtdResidsArray = new double[n][1]; for (int i = 0; i < n; i++) { unwtdResidsArray[i][0] = values[i] - meanVal; } RealMatrix unwtdResids = new BlockRealMatrix(unwtdResidsArray); RealMatrix transUnwtdResids = unwtdResids.transpose(); RealMatrix product = transUnwtdResids.multiply(omega); RealMatrix sumWtdResids = product.multiply(unwtdResids); double mswd = 0.0; double prob = 0.0; if (n > 1) { mswd = sumWtdResids.getEntry(0, 0) / (n - 1); // http://commons.apache.org/proper/commons-math/apidocs/org/apache/commons/math3/distribution/FDistribution.html FDistribution fdist = new org.apache.commons.math3.distribution.FDistribution((n - 1), 1E9); prob = 1.0 - fdist.cumulativeProbability(mswd); } results.setBad(false); results.setMeanVal(meanVal); results.setSigmaMeanVal(meanValSigma); results.setMswd(mswd); results.setProb(prob); } return results; }
From source file:org.easotope.shared.math.ExponentialDecayFitter.java
private void calculateIfNecessary() { if (!needsCalculation) { return;/*from w w w . ja v a 2s. c o m*/ } a = b = c = Double.NaN; if (pointList.size() < 2) { return; } Point[] point = pointList.toArray(new Point[pointList.size()]); Arrays.sort(point, new PointComparator()); double[] s = new double[point.length]; s[0] = 0.0; for (int k = 1; k < point.length; k++) { s[k] = s[k - 1] + 0.5 * (point[k].y + point[k - 1].y) * (point[k].x - point[k - 1].x); } // calculate m RealMatrix m = MatrixUtils.createRealMatrix(2, 2); for (int k = 0; k < point.length; k++) { m.setEntry(0, 0, m.getEntry(0, 0) + (point[k].x - point[0].x) * (point[k].x - point[0].x)); m.setEntry(0, 1, m.getEntry(0, 1) + (point[k].x - point[0].x) * s[k]); m.setEntry(1, 1, m.getEntry(1, 1) + s[k] * s[k]); } m.setEntry(1, 0, m.getEntry(0, 1)); // invert m try { m = MatrixUtils.inverse(m); } catch (Exception e) { return; } // calculate n RealMatrix n = MatrixUtils.createRealMatrix(2, 1); for (int k = 0; k < point.length; k++) { n.setEntry(0, 0, n.getEntry(0, 0) + (point[k].y - point[0].y) * (point[k].x - point[0].x)); n.setEntry(1, 0, n.getEntry(1, 0) + (point[k].y - point[0].y) * s[k]); } // calculate c c = m.multiply(n).getEntry(1, 0); // calculate m m = MatrixUtils.createRealMatrix(2, 2); m.setEntry(0, 0, point.length); for (int k = 0; k < point.length; k++) { m.setEntry(0, 1, m.getEntry(0, 1) + Math.exp(c * point[k].x)); m.setEntry(1, 1, m.getEntry(1, 1) + Math.exp(2.0 * c * point[k].x)); } m.setEntry(1, 0, m.getEntry(0, 1)); // invert m try { m = MatrixUtils.inverse(m); } catch (Exception e) { return; } // calculate n n = MatrixUtils.createRealMatrix(2, 1); for (int k = 0; k < point.length; k++) { n.setEntry(0, 0, n.getEntry(0, 0) + point[k].y); n.setEntry(1, 0, n.getEntry(1, 0) + point[k].y * Math.exp(c * point[k].x)); } // calculate a & b RealMatrix r = m.multiply(n); a = r.getEntry(0, 0); b = r.getEntry(1, 0); }