List of usage examples for org.apache.commons.math3.linear RealVector getEntry
public abstract double getEntry(int index) throws OutOfRangeException;
From source file:edu.stanford.cfuller.imageanalysistools.fitting.GaussianImageObject.java
/** * Fits this object to a 3-dimensional gaussian, and estimates error and goodness of fit. * @param p The parameters for the current analysis. *///from w ww .j ava2s . c om public void fitPosition(ParameterDictionary p) { if (this.sizeInPixels == 0) { this.nullifyImages(); return; } this.fitParametersByChannel = new java.util.ArrayList<FitParameters>(); this.fitR2ByChannel = new java.util.ArrayList<Double>(); this.fitErrorByChannel = new java.util.ArrayList<Double>(); this.nPhotonsByChannel = new java.util.ArrayList<Double>(); GaussianFitter3D gf = new GaussianFitter3D(); //System.out.println(this.parent.getDimensionSizes().getZ()); int numChannels = 0; if (p.hasKey(NUM_WAVELENGTHS_PARAM)) { numChannels = p.getIntValueForKey(NUM_WAVELENGTHS_PARAM); } else { numChannels = this.parent.getDimensionSizes().get(ImageCoordinate.C); } for (int channelIndex = 0; channelIndex < numChannels; channelIndex++) { RealVector fitParameters = new ArrayRealVector(7, 0.0); double ppg = p.getDoubleValueForKey(PHOTONS_PER_LEVEL_PARAM); this.parentBoxMin.set(ImageCoordinate.C, channelIndex); this.parentBoxMax.set(ImageCoordinate.C, channelIndex + 1); this.boxImages(); List<Double> x = new java.util.ArrayList<Double>(); List<Double> y = new java.util.ArrayList<Double>(); List<Double> z = new java.util.ArrayList<Double>(); List<Double> f = new java.util.ArrayList<Double>(); for (ImageCoordinate ic : this.parent) { x.add((double) ic.get(ImageCoordinate.X)); y.add((double) ic.get(ImageCoordinate.Y)); z.add((double) ic.get(ImageCoordinate.Z)); f.add((double) parent.getValue(ic)); } xValues = new double[x.size()]; yValues = new double[y.size()]; zValues = new double[z.size()]; functionValues = new double[f.size()]; double xCentroid = 0; double yCentroid = 0; double zCentroid = 0; double totalCounts = 0; for (int i = 0; i < x.size(); i++) { xValues[i] = x.get(i); yValues[i] = y.get(i); zValues[i] = z.get(i); functionValues[i] = f.get(i) * ppg; xCentroid += xValues[i] * functionValues[i]; yCentroid += yValues[i] * functionValues[i]; zCentroid += zValues[i] * functionValues[i]; totalCounts += functionValues[i]; } xCentroid /= totalCounts; yCentroid /= totalCounts; zCentroid /= totalCounts; //z sometimes seems to be a bit off... trying (20110415) to go back to max value pixel at x,y centroid int xRound = (int) Math.round(xCentroid); int yRound = (int) Math.round(yCentroid); double maxVal = 0; int maxInd = 0; double minZ = Double.MAX_VALUE; double maxZ = 0; for (int i = 0; i < x.size(); i++) { if (zValues[i] < minZ) minZ = zValues[i]; if (zValues[i] > maxZ) maxZ = zValues[i]; if (xValues[i] == xRound && yValues[i] == yRound) { if (functionValues[i] > maxVal) { maxVal = functionValues[i]; maxInd = (int) zValues[i]; } } } zCentroid = maxInd; //parameter ordering: amplitude, var x-y, var z, x/y/z coords, background //amplitude: find the max value; background: find the min value double maxValue = 0; double minValue = Double.MAX_VALUE; for (ImageCoordinate ic : this.parent) { if (parent.getValue(ic) > maxValue) maxValue = parent.getValue(ic); if (parent.getValue(ic) < minValue) minValue = parent.getValue(ic); } fitParameters.setEntry(0, (maxValue - minValue) * 0.95); fitParameters.setEntry(6, minValue + 0.05 * (maxValue - minValue)); //positions fitParameters.setEntry(3, xCentroid); fitParameters.setEntry(4, yCentroid); fitParameters.setEntry(5, zCentroid); //variances final double limitedWidthxy = 200; final double limitedWidthz = 500; double sizex = limitedWidthxy / p.getDoubleValueForKey(PIXELSIZE_PARAM); double sizez = limitedWidthz / p.getDoubleValueForKey(SECTIONSIZE_PARAM); if (p.hasKey(Z_WIDTH_PARAM)) { sizez = p.getDoubleValueForKey(Z_WIDTH_PARAM); } if (p.hasKey(XY_WIDTH_PARAM)) { sizex = p.getDoubleValueForKey(XY_WIDTH_PARAM); } fitParameters.setEntry(1, sizex / 2); fitParameters.setEntry(2, sizez / 2); //amplitude and background are in arbitrary intensity units; convert to photon counts fitParameters.setEntry(0, fitParameters.getEntry(0) * ppg); fitParameters.setEntry(6, fitParameters.getEntry(6) * ppg); //System.out.println("guess: " + fitParameters); //do the fit fitParameters = gf.fit(this, fitParameters, ppg); //System.out.println("fit: " + fitParameters); FitParameters fp = new FitParameters(); fp.setPosition(ImageCoordinate.X, fitParameters.getEntry(3)); fp.setPosition(ImageCoordinate.Y, fitParameters.getEntry(4)); fp.setPosition(ImageCoordinate.Z, fitParameters.getEntry(5)); fp.setSize(ImageCoordinate.X, fitParameters.getEntry(1)); fp.setSize(ImageCoordinate.Y, fitParameters.getEntry(1)); fp.setSize(ImageCoordinate.Z, fitParameters.getEntry(2)); fp.setAmplitude(fitParameters.getEntry(0)); fp.setBackground(fitParameters.getEntry(6)); fitParametersByChannel.add(fp); //calculate R2 double residualSumSquared = 0; double mean = 0; double variance = 0; double R2 = 0; double n_photons = 0; for (int i = 0; i < this.xValues.length; i++) { residualSumSquared += Math.pow(GaussianFitter3D.fitResidual(functionValues[i], xValues[i], yValues[i], zValues[i], fitParameters), 2); mean += functionValues[i]; n_photons += functionValues[i] - fitParameters.getEntry(6); } mean /= functionValues.length; for (int i = 0; i < this.xValues.length; i++) { variance += Math.pow(functionValues[i] - mean, 2); } R2 = 1 - (residualSumSquared / variance); this.fitR2ByChannel.add(R2); this.unboxImages(); //calculate fit error double s_xy = fitParameters.getEntry(1) * fitParameters.getEntry(1) * Math.pow(p.getDoubleValueForKey(PIXELSIZE_PARAM), 2); double s_z = fitParameters.getEntry(2) * fitParameters.getEntry(2) * Math.pow(p.getDoubleValueForKey(SECTIONSIZE_PARAM), 2); //s_z = 0; //remove!! double error = (2 * s_xy + s_z) / (n_photons - 1);// + 4*Math.sqrt(Math.PI) * Math.pow(2*s_xy,1.5)*Math.pow(fitParameters.getEntry(6),2)/(p.getDoubleValueForKey("pixelsize_nm")*n_photons*n_photons); double b = fitParameters.getEntry(6); double a = p.getDoubleValueForKey(PIXELSIZE_PARAM); double alpha = p.getDoubleValueForKey(SECTIONSIZE_PARAM); double sa_x = s_xy + Math.pow(a, 2) / 12; double sa_z = s_z + Math.pow(alpha, 2) / 12; //System.out.printf("b = %f, a = %f, alpha = %f, s_xy = %f, s_z = %f, n= %f\n", b, a, alpha, s_xy, s_z, n_photons); double error_x = sa_x / n_photons * (16.0 / 9.0 + 8 * Math.PI * sa_x * b * b / (n_photons * Math.pow(p.getDoubleValueForKey(PIXELSIZE_PARAM), 2))); double error_z = sa_z / n_photons * (16.0 / 9.0 + 8 * Math.PI * sa_z * b * b / (n_photons * Math.pow(p.getDoubleValueForKey(SECTIONSIZE_PARAM), 2))); double A = 1.0 / (2 * Math.sqrt(2) * Math.pow(Math.PI, 1.5) * Math.sqrt(sa_z) * sa_x); ErrIntFunc eif = new ErrIntFunc(); eif.setParams(b, n_photons, A, sa_z, sa_x, a, alpha); LegendreGaussIntegrator lgi = new LegendreGaussIntegrator(5, 10, 1000); //integrate over 10*width of PSF in z double size = 10 * Math.sqrt(sa_z); double intpart = 0; try { if (b < 0) throw new ConvergenceException(new DummyLocalizable("negative background!")); // a negative value for b seems to cause the integration to hang, preventing the program from progressing intpart = lgi.integrate(10000, eif, -size, size); double fullIntPart = intpart + Math.pow(2 * Math.PI, 1.5) * sa_x * A / Math.sqrt(sa_z); error_x = Math.sqrt(2 / (n_photons * sa_z / (2 * sa_z + sa_x) * fullIntPart)); error_z = Math.sqrt(2 / (n_photons * sa_x / (2 * sa_z + sa_x) * fullIntPart)); } catch (ConvergenceException e) { LoggingUtilities.getLogger().severe("Integration error: " + e.getMessage()); error_x = -1; error_z = -1; } if (error_x > 0 && error_z > 0) { error = Math.sqrt(2 * error_x * error_x + error_z * error_z); } else { error = Double.NaN; } this.fitErrorByChannel.add(error); this.positionsByChannel.add(fitParameters.getSubVector(3, 3)); this.nPhotonsByChannel.add(n_photons); } this.hadFittingError = false; }
From source file:edu.stanford.cfuller.imageanalysistools.fitting.GaussianImageObjectWithCovariance.java
/** * Fits this object to a 3-dimensional gaussian, and estimates error and goodness of fit. * @param p The parameters for the current analysis. *///from w w w . j a v a2 s. c om public void fitPosition(ParameterDictionary p) { if (this.sizeInPixels == 0) { this.nullifyImages(); return; } this.fitParametersByChannel = new java.util.ArrayList<FitParameters>(); this.fitR2ByChannel = new java.util.ArrayList<Double>(); this.fitErrorByChannel = new java.util.ArrayList<Double>(); this.nPhotonsByChannel = new java.util.ArrayList<Double>(); GaussianFitter3DWithCovariance gf = new GaussianFitter3DWithCovariance(); //System.out.println(this.parent.getDimensionSizes().getZ()); int numChannels = 0; if (p.hasKey("num_wavelengths")) { numChannels = p.getIntValueForKey("num_wavelengths"); } else { numChannels = this.parent.getDimensionSizes().get(ImageCoordinate.C); } //for (int channelIndex = 0; channelIndex < this.parent.getDimensionSizes().getC(); channelIndex++) { for (int channelIndex = 0; channelIndex < numChannels; channelIndex++) { RealVector fitParameters = new ArrayRealVector(11, 0.0); double ppg = p.getDoubleValueForKey("photons_per_greylevel"); this.parentBoxMin.set(ImageCoordinate.C, channelIndex); this.parentBoxMax.set(ImageCoordinate.C, channelIndex + 1); this.boxImages(); List<Double> x = new java.util.ArrayList<Double>(); List<Double> y = new java.util.ArrayList<Double>(); List<Double> z = new java.util.ArrayList<Double>(); List<Double> f = new java.util.ArrayList<Double>(); for (ImageCoordinate ic : this.parent) { x.add((double) ic.get(ImageCoordinate.X)); y.add((double) ic.get(ImageCoordinate.Y)); z.add((double) ic.get(ImageCoordinate.Z)); f.add((double) parent.getValue(ic)); } xValues = new double[x.size()]; yValues = new double[y.size()]; zValues = new double[z.size()]; functionValues = new double[f.size()]; double xCentroid = 0; double yCentroid = 0; double zCentroid = 0; double totalCounts = 0; for (int i = 0; i < x.size(); i++) { xValues[i] = x.get(i); yValues[i] = y.get(i); zValues[i] = z.get(i); functionValues[i] = f.get(i) * ppg; xCentroid += xValues[i] * functionValues[i]; yCentroid += yValues[i] * functionValues[i]; zCentroid += zValues[i] * functionValues[i]; totalCounts += functionValues[i]; } xCentroid /= totalCounts; yCentroid /= totalCounts; zCentroid /= totalCounts; //z sometimes seems to be a bit off... trying (20110415) to go back to max value pixel at x,y centroid int xRound = (int) Math.round(xCentroid); int yRound = (int) Math.round(yCentroid); double maxVal = 0; int maxInd = 0; double minZ = Double.MAX_VALUE; double maxZ = 0; for (int i = 0; i < x.size(); i++) { if (zValues[i] < minZ) minZ = zValues[i]; if (zValues[i] > maxZ) maxZ = zValues[i]; if (xValues[i] == xRound && yValues[i] == yRound) { if (functionValues[i] > maxVal) { maxVal = functionValues[i]; maxInd = (int) zValues[i]; } } } zCentroid = maxInd; //parameter ordering: amplitude, var x-y, var z, x/y/z coords, background //amplitude: find the max value; background: find the min value double maxValue = 0; double minValue = Double.MAX_VALUE; for (ImageCoordinate ic : this.parent) { if (parent.getValue(ic) > maxValue) maxValue = parent.getValue(ic); if (parent.getValue(ic) < minValue) minValue = parent.getValue(ic); } fitParameters.setEntry(0, (maxValue - minValue) * 0.95); fitParameters.setEntry(10, minValue + 0.05 * (maxValue - minValue)); //positions fitParameters.setEntry(7, xCentroid); fitParameters.setEntry(8, yCentroid); fitParameters.setEntry(9, zCentroid); //variances final double limitedWidthxy = 200; final double limitedWidthz = 500; double sizex = limitedWidthxy / p.getDoubleValueForKey("pixelsize_nm"); double sizey = limitedWidthxy / p.getDoubleValueForKey("pixelsize_nm"); double sizez = limitedWidthz / p.getDoubleValueForKey("z_sectionsize_nm"); if (p.hasKey(Z_WIDTH_PARAM)) { sizez = p.getDoubleValueForKey(Z_WIDTH_PARAM); } if (p.hasKey(XY_WIDTH_PARAM)) { sizex = p.getDoubleValueForKey(XY_WIDTH_PARAM); sizey = p.getDoubleValueForKey(XY_WIDTH_PARAM); } fitParameters.setEntry(1, sizex / 2); fitParameters.setEntry(2, sizey / 2); fitParameters.setEntry(3, sizez / 2); //covariances -- guess zero to start fitParameters.setEntry(4, 0.0); fitParameters.setEntry(5, 0.0); fitParameters.setEntry(6, 0.0); //amplitude and background are in arbitrary intensity units; convert to photon counts fitParameters.setEntry(0, fitParameters.getEntry(0) * ppg); fitParameters.setEntry(10, fitParameters.getEntry(10) * ppg); //System.out.println("guess: " + fitParameters); //do the fit //System.out.println("Initial for object " + this.label + ": " + fitParameters.toString()); fitParameters = gf.fit(this, fitParameters, ppg); //System.out.println("Parameters for object " + this.label + ": " + fitParameters.toString()); //System.out.println("fit: " + fitParameters); FitParameters fp = new FitParameters(); fp.setPosition(ImageCoordinate.X, fitParameters.getEntry(7)); fp.setPosition(ImageCoordinate.Y, fitParameters.getEntry(8)); fp.setPosition(ImageCoordinate.Z, fitParameters.getEntry(9)); fp.setSize(ImageCoordinate.X, fitParameters.getEntry(1)); fp.setSize(ImageCoordinate.Y, fitParameters.getEntry(2)); fp.setSize(ImageCoordinate.Z, fitParameters.getEntry(3)); fp.setAmplitude(fitParameters.getEntry(0)); fp.setBackground(fitParameters.getEntry(10)); fp.setOtherParameter("corr_xy", fitParameters.getEntry(4)); fp.setOtherParameter("corr_xz", fitParameters.getEntry(5)); fp.setOtherParameter("corr_yz", fitParameters.getEntry(6)); fitParametersByChannel.add(fp); // System.out.println(fitParameters); //calculate R2 double residualSumSquared = 0; double mean = 0; double variance = 0; double R2 = 0; double n_photons = 0; for (int i = 0; i < this.xValues.length; i++) { residualSumSquared += Math.pow(GaussianFitter3DWithCovariance.fitResidual(functionValues[i], xValues[i], yValues[i], zValues[i], fitParameters), 2); mean += functionValues[i]; n_photons += functionValues[i] - fitParameters.getEntry(10); } mean /= functionValues.length; for (int i = 0; i < this.xValues.length; i++) { variance += Math.pow(functionValues[i] - mean, 2); } R2 = 1 - (residualSumSquared / variance); this.fitR2ByChannel.add(R2); this.unboxImages(); //calculate fit error -- these are wrong for the case with unequal x-y variances and covariance allowed, but leave them for now. double s_xy = fitParameters.getEntry(1) * fitParameters.getEntry(1) * Math.pow(p.getDoubleValueForKey("pixelsize_nm"), 2); double s_z = fitParameters.getEntry(3) * fitParameters.getEntry(3) * Math.pow(p.getDoubleValueForKey("z_sectionsize_nm"), 2); //s_z = 0; //remove!! double error = (2 * s_xy + s_z) / (n_photons - 1);// + 4*Math.sqrt(Math.PI) * Math.pow(2*s_xy,1.5)*Math.pow(fitParameters.getEntry(6),2)/(p.getDoubleValueForKey("pixelsize_nm")*n_photons*n_photons); double b = fitParameters.getEntry(10); double a = p.getDoubleValueForKey("pixelsize_nm"); double alpha = p.getDoubleValueForKey("z_sectionsize_nm"); double sa_x = s_xy + Math.pow(a, 2) / 12; double sa_z = s_z + Math.pow(alpha, 2) / 12; double error_x = sa_x / n_photons * (16.0 / 9.0 + 8 * Math.PI * sa_x * b * b / (n_photons * Math.pow(p.getDoubleValueForKey("pixelsize_nm"), 2))); double error_z = sa_z / n_photons * (16.0 / 9.0 + 8 * Math.PI * sa_z * b * b / (n_photons * Math.pow(p.getDoubleValueForKey("z_sectionsize_nm"), 2))); double A = 1.0 / (2 * Math.sqrt(2) * Math.pow(Math.PI, 1.5) * Math.sqrt(sa_z) * sa_x); ErrIntFunc eif = new ErrIntFunc(); eif.setParams(b, n_photons, A, sa_z, sa_x, a, alpha); //LegendreGaussIntegrator lgi = new LegendreGaussIntegrator(5, 10, 1000); //integrate over 10*width of PSF in z //double size = 10*Math.sqrt(sa_z); //double intpart = 0; // try { // // if (b < 0) throw new ConvergenceException(new DummyLocalizable("negative background!")); // a negative value for b seems to cause the integration to hang, preventing the program from progressing // // intpart = lgi.integrate(10000, eif, -size, size); // // double fullIntPart = intpart + Math.pow(2*Math.PI, 1.5)*sa_x*A/Math.sqrt(sa_z); // // error_x = Math.sqrt(2/(n_photons*sa_z/(2*sa_z + sa_x)*fullIntPart)); // error_z = Math.sqrt(2/(n_photons*sa_x/(2*sa_z + sa_x)*fullIntPart)); // // } catch (ConvergenceException e) { // LoggingUtilities.getLogger().severe("Integration error: " + e.getMessage()); // error_x = -1; // error_z = -1; // } // if (error_x > 0 && error_z > 0) { error = Math.sqrt(2 * error_x * error_x + error_z * error_z); } else { error = Double.NaN; } error = 0; // until a better error calculation this.fitErrorByChannel.add(error); this.positionsByChannel.add(fitParameters.getSubVector(7, 3)); this.nPhotonsByChannel.add(n_photons); } this.hadFittingError = false; this.nullifyImages(); }
From source file:com.joptimizer.optimizers.LPPresolver.java
/** * This method is just for testing scope. *//*from ww w. j av a 2 s . c om*/ private void checkProgress(DoubleMatrix1D c, DoubleMatrix2D A, DoubleMatrix1D b, DoubleMatrix1D lb, DoubleMatrix1D ub, DoubleMatrix1D ylb, DoubleMatrix1D yub, DoubleMatrix1D zlb, DoubleMatrix1D zub) { if (this.expectedSolution == null) { return; } if (Double.isNaN(this.expectedTolerance)) { //for this to work properly, this method must be called at least one time before presolving operations start RealVector X = MatrixUtils.createRealVector(expectedSolution); RealMatrix AMatrix = MatrixUtils.createRealMatrix(A.toArray()); RealVector Bvector = MatrixUtils.createRealVector(b.toArray()); RealVector Axb = AMatrix.operate(X).subtract(Bvector); double norm = Axb.getNorm(); this.expectedTolerance = Math.max(1.e-7, 1.5 * norm); } double tolerance = this.expectedTolerance; log.debug("tolerance: " + tolerance); RealVector X = MatrixUtils.createRealVector(expectedSolution); RealMatrix AMatrix = MatrixUtils.createRealMatrix(A.toArray()); RealVector Bvector = MatrixUtils.createRealVector(b.toArray()); //logger.debug("A.X-b: " + ArrayUtils.toString(originalA.operate(X).subtract(originalB))); //nz rows for (int i = 0; i < vRowPositions.length; i++) { short[] vRowPositionsI = vRowPositions[i]; for (short nzJ : vRowPositionsI) { if (Double.compare(A.getQuick(i, nzJ), 0.) == 0) { log.debug("entry " + i + "," + nzJ + " est zero: " + A.getQuick(i, nzJ)); throw new IllegalStateException(); } } } //nz columns for (int j = 0; j < vColPositions.length; j++) { short[] vColPositionsJ = vColPositions[j]; for (short nzI : vColPositionsJ) { if (Double.compare(A.getQuick(nzI, j), 0.) == 0) { log.debug("entry (" + nzI + "," + j + ") est zero: " + A.getQuick(nzI, j)); throw new IllegalStateException(); } } } //nz Aij for (int i = 0; i < A.rows(); i++) { short[] vRowPositionsI = vRowPositions[i]; for (int j = 0; j < A.columns(); j++) { if (Double.compare(Math.abs(A.getQuick(i, j)), 0.) != 0) { if (!ArrayUtils.contains(vRowPositionsI, (short) j)) { log.debug("entry " + i + "," + j + " est non-zero: " + A.getQuick(i, j)); throw new IllegalStateException(); } if (!ArrayUtils.contains(vColPositions[j], (short) i)) { log.debug("entry " + i + "," + j + " est non-zero: " + A.getQuick(i, j)); throw new IllegalStateException(); } } } } // //boolean deepCheckA = true; // boolean deepCheckA = false; // if(deepCheckA){ // //check for 0-rows // List<Integer> zeroRows = new ArrayList<Integer>(); // for(int i=0; i<A.rows(); i++){ // boolean isNotZero = false; // for(int j=0;!isNotZero && j<A.columns(); j++){ // isNotZero = Double.compare(0., A.getQuick(i, j))!=0; // } // if(!isNotZero){ // zeroRows.add(zeroRows.size(), i); // } // } // if(!zeroRows.isEmpty()){ // log.debug("All 0 entries in rows " + ArrayUtils.toString(zeroRows)); // //log.debug(ArrayUtils.toString(A.toArray())); // throw new IllegalStateException(); // } // // //check for 0-columns // List<Integer> zeroCols = new ArrayList<Integer>(); // for(int j=0; j<A.columns(); j++){ // boolean isNotZero = false; // for(int i=0;!isNotZero && i<A.rows(); i++){ // isNotZero = Double.compare(0., A.getQuick(i, j))!=0; // } // if(!isNotZero){ // zeroCols.add(zeroCols.size(), j); // } // } // if(!zeroCols.isEmpty()){ // log.debug("All 0 entries in columns " + ArrayUtils.toString(zeroCols)); // //log.debug(ArrayUtils.toString(A.toArray())); // throw new IllegalStateException(); // } // // // check rank(A): must be A pXn with rank(A)=p < n // QRSparseFactorization qr = null; // boolean factOK = true; // try{ // qr = new QRSparseFactorization((SparseDoubleMatrix2D)A); // qr.factorize(); // }catch(Exception e){ // factOK = false; // log.warn("Warning", e); // } // if(factOK){ // log.debug("p : " + AMatrix.getRowDimension()); // log.debug("n : " + AMatrix.getColumnDimension()); // log.debug("full rank: " + qr.hasFullRank()); // if(!(A.rows() < A.columns())){ // log.debug("!( p < n )"); // throw new IllegalStateException(); // } // if(!qr.hasFullRank()){ // log.debug("not full rank A matrix"); // throw new IllegalStateException(); // } // } // } //A.x = b RealVector Axb = AMatrix.operate(X).subtract(Bvector); double norm = Axb.getNorm(); log.debug("|| A.x-b ||: " + norm); if (norm > tolerance) { //where is the error? for (int i = 0; i < Axb.getDimension(); i++) { if (Math.abs(Axb.getEntry(i)) > tolerance) { log.debug("entry " + i + ": " + Axb.getEntry(i)); throw new IllegalStateException(); } } throw new IllegalStateException(); } //upper e lower for (int i = 0; i < X.getDimension(); i++) { if (X.getEntry(i) + tolerance < lb.getQuick(i)) { log.debug("lower bound " + i + " not respected: lb=" + lb.getQuick(i) + ", value=" + X.getEntry(i)); throw new IllegalStateException(); } if (X.getEntry(i) > ub.getQuick(i) + tolerance) { log.debug("upper bound " + i + " not respected: ub=" + ub.getQuick(i) + ", value=" + X.getEntry(i)); throw new IllegalStateException(); } } }
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 ww w.j a v a 2 s.c o m // 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.apache.predictionio.examples.java.recommendations.tutorial1.Algorithm.java
@Override public Float predict(Model model, Query query) { RealVector itemVector = model.itemSimilarity.get(query.iid); RealVector userVector = model.userHistory.get(query.uid); if (itemVector == null) { // cold start item, can't be handled by this algo, return hard code value. return Float.NaN; } else if (userVector == null) { // new user, can't be handled by this algo, return hard code value. return Float.NaN; } else {/*from ww w . j a v a 2 s . c o m*/ //logger.info("(" + query.uid + "," + query.iid + ")"); //logger.info(itemVector.toString()); //logger.info(userVector.toString()); double accum = 0.0; double accumSim = 0.0; for (int i = 0; i < itemVector.getDimension(); i++) { double weight = itemVector.getEntry(i); double rating = userVector.getEntry(i); if ((weight != 0) && (rating != 0)) { accum += weight * rating; accumSim += Math.abs(weight); } } if (accumSim == 0.0) { return Float.NaN; } else { return (float) (accum / accumSim); } } }
From source file:org.apache.predictionio.examples.java.recommendations.tutorial4.CollaborativeFilteringAlgorithm.java
@Override public Float predict(CollaborativeFilteringModel model, Query query) { RealVector itemVector = model.itemSimilarity.get(query.iid); RealVector userVector = model.userHistory.get(query.uid); if (itemVector == null) { // cold start item, can't be handled by this algo, return hard code value. return Float.NaN; } else if (userVector == null) { // new user, can't be handled by this algo, return hard code value. return Float.NaN; } else {//from w w w .jav a 2 s .c om //logger.info("(" + query.uid + "," + query.iid + ")"); //logger.info(itemVector.toString()); //logger.info(userVector.toString()); double accum = 0.0; double accumSim = 0.0; for (int i = 0; i < itemVector.getDimension(); i++) { double weight = itemVector.getEntry(i); double rating = userVector.getEntry(i); if ((weight != 0) && (rating != 0)) { accum += weight * rating; accumSim += Math.abs(weight); } } if (accumSim == 0.0) { return Float.NaN; } else { return (float) (accum / accumSim); } } }
From source file:org.briljantframework.array.Matrices.java
/** * View the double array as a {@link RealVector}. * * @param array the array/* w ww. j av a 2s.c om*/ * @return a real vector view */ public static RealVector asRealVector(DoubleArray array) { Check.argument(array.isVector(), CAN_ONLY_VIEW_1D_ARRAYS); return new RealVector() { @Override public int getDimension() { return array.size(); } @Override public double getEntry(int index) throws OutOfRangeException { return array.get(index); } @Override public void setEntry(int index, double value) throws OutOfRangeException { array.set(index, value); } @Override public RealVector append(RealVector v) { ArrayRealVector vector = new ArrayRealVector(v.getDimension() + array.size()); copyFromArray(vector); for (int i = 0; i < v.getDimension(); i++) { vector.setEntry(i, v.getEntry(i)); } return vector; } private void copyFromArray(ArrayRealVector vector) { for (int i = 0; i < array.size(); i++) { vector.setEntry(i, array.get(i)); } } @Override public RealVector append(double d) { ArrayRealVector vector = new ArrayRealVector(array.size() + 1); copyFromArray(vector); vector.setEntry(array.size(), d); return vector; } @Override public RealVector getSubVector(int index, int n) throws NotPositiveException, OutOfRangeException { return asRealVector(array.get(Range.of(index, index + n))); } @Override public void setSubVector(int index, RealVector v) throws OutOfRangeException { for (int i = 0; i < array.size(); i++) { array.set(i + index, v.getEntry(i)); } } @Override public boolean isNaN() { return Arrays.any(array, Double::isNaN); } @Override public boolean isInfinite() { return Arrays.any(array, Double::isInfinite); } @Override public RealVector copy() { return asRealVector(array.copy()); } @Override @Deprecated public RealVector ebeDivide(RealVector v) throws DimensionMismatchException { if (v.getDimension() != array.size()) { throw new DimensionMismatchException(v.getDimension(), array.size()); } RealVector a = new ArrayRealVector(array.size()); for (int i = 0; i < array.size(); i++) { a.setEntry(i, array.get(i) / v.getEntry(i)); } return a; } @Override @Deprecated public RealVector ebeMultiply(RealVector v) throws DimensionMismatchException { if (v.getDimension() != array.size()) { throw new DimensionMismatchException(v.getDimension(), array.size()); } RealVector a = new ArrayRealVector(array.size()); for (int i = 0; i < array.size(); i++) { a.setEntry(i, array.get(i) * v.getEntry(i)); } return a; } }; }
From source file:org.eclipse.dataset.LinearAlgebra.java
private static Dataset createDataset(RealVector v) { DoubleDataset r = new DoubleDataset(v.getDimension()); int size = r.getSize(); if (v instanceof ArrayRealVector) { double[] data = ((ArrayRealVector) v).getDataRef(); for (int i = 0; i < size; i++) { r.setAbs(i, data[i]);//from ww w .j a v a 2s . c o m } } else { for (int i = 0; i < size; i++) { r.setAbs(i, v.getEntry(i)); } } return r; }
From source file:org.eclipse.january.dataset.LinearAlgebra.java
private static Dataset createDataset(RealVector v) { DoubleDataset r = DatasetFactory.zeros(DoubleDataset.class, v.getDimension()); int size = r.getSize(); if (v instanceof ArrayRealVector) { double[] data = ((ArrayRealVector) v).getDataRef(); for (int i = 0; i < size; i++) { r.setAbs(i, data[i]);/* ww w . j a va2s. c om*/ } } else { for (int i = 0; i < size; i++) { r.setAbs(i, v.getEntry(i)); } } return r; }
From source file:org.grouplens.lenskit.diffusion.Iterative.IterativeDiffusionItemScorer.java
/** * Score items by computing predicted ratings. * * @see ItemScoreAlgorithm#scoreItems(ItemItemModel, org.grouplens.lenskit.vectors.SparseVector, org.grouplens.lenskit.vectors.MutableSparseVector, NeighborhoodScorer) *//*from ww w . ja v a 2 s. c o m*/ @Override public void score(long user, @Nonnull MutableSparseVector scores) { UserHistory<? extends Event> history = dao.getEventsForUser(user, summarizer.eventTypeWanted()); if (history == null) { history = History.forUser(user); } SparseVector summary = summarizer.summarize(history); VectorTransformation transform = normalizer.makeTransformation(user, summary); MutableSparseVector normed = summary.mutableCopy(); transform.apply(normed); scores.clear(); int numItems = 1682; //algorithm.scoreItems(model, normed, scores, scorer); int num_updates = 300; double update_rate = 1; double threshold = 0.01; RealVector z_out = diffusionMatrix.preMultiply(VectorUtils.toRealVector(numItems, normed)); boolean updated = true; LongSortedSet known = normed.keySet(); int count_iter = 0; for (int i = 0; i < num_updates && updated; i++) { updated = false; RealVector temp = diffusionMatrix.preMultiply(z_out); temp.mapMultiplyToSelf(z_out.getNorm() / temp.getNorm()); RealVector temp_diff = z_out.add(temp.mapMultiplyToSelf(-1.0)); for (int j = 0; j < numItems; j++) { if (!known.contains((long) (j + 1))) { //if the rating is not one of the known ones if (Math.abs(temp_diff.getEntry(j)) > threshold) { // if difference is large enough, update updated = true; z_out.setEntry(j, (1.0 - update_rate) * z_out.getEntry(j) + update_rate * temp.getEntry(j)); } } } count_iter++; } System.out.println(count_iter); LongSortedSet testDomain = scores.keyDomain(); //fill up the score vector for (int i = 0; i < numItems; i++) { if (testDomain.contains((long) (i + 1))) { scores.set((long) (i + 1), z_out.getEntry(i)); } } // untransform the scores transform.unapply(scores); System.out.println(scores); }