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.NelderMeadMinimizer.java
/** * Constructs the initial simplex that is the starting point of the optimization given an initial guess at the minimum. * * This method will attempt to guess the scale of each parameters, but it is preferable to specify this manually using the other form of * generateInitialSimplex if any information about these scales is known. * * @param initialPoint The initial guess at the minimum, one component per parameter. * @return A matrix containing p + 1 rows, each of which is one set of p parameters, which specify the simplex. *///from w ww . j ava 2s . co m public RealMatrix generateInitialSimplex(RealVector initialPoint) { final double constantScale = 0.1; RealVector componentScales = initialPoint.mapMultiply(constantScale); //if the initial point has zeros in it, those entries will not be optimized //perturb slightly to allow optimization for (int i = 0; i < componentScales.getDimension(); i++) { if (componentScales.getEntry(i) == 0.0) { componentScales.setEntry(i, constantScale); } } return generateInitialSimplex(initialPoint, componentScales); }
From source file:edu.stanford.cfuller.imageanalysistools.filter.LocalMaximumSeparabilityThresholdingFilter.java
/** * Applies the filter to an Image, optionally turning on adaptive determination of the increment size used to find the threshold, and specifying a size for the threshold determination increment. * Turning on adaptive determination of the increment will generally make the threshold slightly less optimal, but can sometimes speed up the filtering, especially * for images with a large dynamic range. * <p>/*w ww .ja va2s . c o m*/ * The increment size specified (in greylevels) will be used to determine the threshold only if adaptive determination is off; otherwise this parameter will be ignored. * @param im The Image to be thresholded; this will be overwritten by the thresholded Image. * @param adaptiveincrement true to turn on adaptive determination of the threshold increment; false to turn it off and use the default value * @param increment the increment size (in greylevels) to use for determining the threshold; must be positive. */ public void apply_ext(WritableImage im, boolean adaptiveincrement, int increment) { Histogram h = new Histogram(im); int thresholdValue = 0; final int numSteps = 1000; double best_eta = 0; int best_index = Integer.MAX_VALUE; int nonzerocounts = h.getTotalCounts() - h.getCounts(0); double meannonzero = h.getMeanNonzero(); ArrayRealVector omega_v = new ArrayRealVector(h.getMaxValue()); ArrayRealVector mu_v = new ArrayRealVector(h.getMaxValue()); ArrayRealVector eta_v = new ArrayRealVector(h.getMaxValue()); int c = 0; if (adaptiveincrement) { increment = (int) ((h.getMaxValue() - h.getMinValue() + 1) * 1.0 / numSteps); if (increment < 1) increment = 1; } for (int k = h.getMinValue(); k < h.getMaxValue() + 1; k += increment) { if (k == 0) continue; omega_v.setEntry(c, h.getCumulativeCounts(k) * 1.0 / nonzerocounts); if (c == 0) { mu_v.setEntry(c, k * omega_v.getEntry(c)); } else { mu_v.setEntry(c, mu_v.getEntry(c - 1) + k * h.getCounts(k) * 1.0 / nonzerocounts); for (int i = k - increment + 1; i < k; i++) { mu_v.setEntry(c, mu_v.getEntry(c) + h.getCounts(i) * i * 1.0 / nonzerocounts); } } double omega = omega_v.getEntry(c); double mu = mu_v.getEntry(c); if (omega > 1e-8 && 1 - omega > 1e-8) { double eta = omega * (1 - omega) * Math.pow((meannonzero - mu) / (1 - omega) - mu / omega, 2); eta_v.setEntry(c, eta); if (eta >= best_eta) { best_eta = eta; best_index = k; } } else { eta_v.setEntry(c, 0); } c++; } int orig_method_best_index = best_index; c = 1; ArrayList<Integer> maxima = new ArrayList<Integer>(); Map<Integer, Integer> k_by_c = new HashMap<Integer, Integer>(); Map<Integer, Integer> c_by_k = new HashMap<Integer, Integer>(); for (int k = h.getMinValue() + 1; k < h.getMaxValue(); k += increment) { //detect if this is a local maximum k_by_c.put(c, k); c_by_k.put(k, c); int lastEntryNotEqual = c - 1; int nextEntryNotEqual = c + 1; while (lastEntryNotEqual > 0 && eta_v.getEntry(lastEntryNotEqual) == eta_v.getEntry(c)) { --lastEntryNotEqual; } while (nextEntryNotEqual < (eta_v.getDimension() - 1) && eta_v.getEntry(nextEntryNotEqual) == eta_v.getEntry(c)) { ++nextEntryNotEqual; } if (eta_v.getEntry(c) > eta_v.getEntry(lastEntryNotEqual) && eta_v.getEntry(c) > eta_v.getEntry(nextEntryNotEqual)) { maxima.add(k); } c++; } //now that we have maxima, try doing a gaussian fit to find the positions. If there's only one, we need to guess at a second RealVector parameters = new ArrayRealVector(6, 0.0); int position0 = 0; int position1 = h.getMaxValue(); if (maxima.size() > 1) { double best_max = 0; double second_best_max = 0; int best_pos = 0; int second_best_pos = 0; for (Integer k : maxima) { int ck = c_by_k.get(k); double eta_k = eta_v.getEntry(ck); if (eta_k > best_max) { second_best_max = best_max; second_best_pos = best_pos; best_max = eta_k; best_pos = ck; } else if (eta_k > second_best_max) { second_best_max = eta_k; second_best_pos = ck; } } position0 = best_pos; position1 = second_best_pos; } else { position0 = c_by_k.get(maxima.get(0)); position1 = (eta_v.getDimension() - position0) / 2 + position0; } //make sure that position 1 is larger than position 0 if (position1 < position0) { int temp = position0; position0 = position1; position1 = temp; } double s = (position1 - position0) / 4.0; parameters.setEntry(0, eta_v.getEntry(position0));//*Math.sqrt(2*Math.PI)*s); parameters.setEntry(1, position0); parameters.setEntry(2, s); parameters.setEntry(3, eta_v.getEntry(position1));//*Math.sqrt(2*Math.PI)*s); parameters.setEntry(4, position1); parameters.setEntry(5, s); DoubleGaussianObjectiveFunction dgof = new DoubleGaussianObjectiveFunction(); dgof.setEta(eta_v); NelderMeadMinimizer nmm = new NelderMeadMinimizer(); RealVector result = nmm.optimize(dgof, parameters); best_index = (int) result.getEntry(4); if (k_by_c.containsKey(best_index)) { best_index = k_by_c.get(best_index); } else { //fall back to the normal global maximum if the fitting seems to have found an invalid value. best_index = orig_method_best_index; } thresholdValue = best_index; if (thresholdValue == Integer.MAX_VALUE) { thresholdValue = 0; } for (ImageCoordinate coord : im) { if (im.getValue(coord) < thresholdValue) im.setValue(coord, 0); } }
From source file:edu.duke.cs.osprey.tupexp.IterativeCGTupleFitter.java
double[] calcFitVals(RealVector rv) { double[] fitVals = new double[numSamp]; for (int s = 0; s < numSamp; s++) { //if(isSampleRestrained(s)){//we need these for validation, if not for AtA aplication ArrayList<Integer> sampTup = tupIndMat.calcSampleTuples(samples.get(s)); for (int t : sampTup) fitVals[s] += rv.getEntry(t); //}//from w w w. j av a2s .com } return fitVals; }
From source file:com.datumbox.framework.core.machinelearning.regression.MatrixLinearRegression.java
/** {@inheritDoc} */ @Override//from w w w .ja va 2 s .c o m protected void _predict(Dataframe newData) { //read model params ModelParameters modelParameters = knowledgeBase.getModelParameters(); Map<Object, Double> thitas = modelParameters.getThitas(); Map<Object, Integer> featureIds = modelParameters.getFeatureIds(); int d = thitas.size(); RealVector coefficients = new OpenMapRealVector(d); for (Map.Entry<Object, Double> entry : thitas.entrySet()) { Integer featureId = featureIds.get(entry.getKey()); coefficients.setEntry(featureId, entry.getValue()); } Map<Integer, Integer> recordIdsReference = new HashMap<>(); //use a mapping between recordIds and rowIds in Matrix DataframeMatrix matrixDataset = DataframeMatrix.parseDataset(newData, recordIdsReference, featureIds); RealMatrix X = matrixDataset.getX(); RealVector Y = X.operate(coefficients); for (Map.Entry<Integer, Record> e : newData.entries()) { Integer rId = e.getKey(); Record r = e.getValue(); int rowId = recordIdsReference.get(rId); newData._unsafe_set(rId, new Record(r.getX(), r.getY(), Y.getEntry(rowId), r.getYPredictedProbabilities())); } //recordIdsReference = null; //matrixDataset = null; }
From source file:com.github.thorbenlindhauer.factor.GaussianFactorTest.java
@Test public void testValueObservation() { GaussianFactor reducedVector = abcFactor.observation(newScope(new ContinuousVariable("A")), new double[] { 2.5d }); // then// www.j a v a 2 s . c om Collection<Variable> newVariables = reducedVector.getVariables().getVariables(); assertThat(newVariables).hasSize(2); assertThat(newVariables).contains(new ContinuousVariable("B"), new ContinuousVariable("C")); // B C A //6.0d, 7.0d, 3.0d //3.0d, 5.5d, 10.0d //4.0d, 6.0d, 3.0d // // X = {B, C}, Y = {A} // precision matrix: K_xx RealMatrix precisionMatrix = reducedVector.getPrecisionMatrix(); assertThat(precisionMatrix.isSquare()).isTrue(); assertThat(precisionMatrix.getColumnDimension()).isEqualTo(2); double precision = precisionMatrix.getEntry(0, 0); assertThat(precision).isEqualTo(6.0d, TestConstants.DOUBLE_VALUE_TOLERANCE); precision = precisionMatrix.getEntry(0, 1); assertThat(precision).isEqualTo(7.0d, TestConstants.DOUBLE_VALUE_TOLERANCE); precision = precisionMatrix.getEntry(1, 0); assertThat(precision).isEqualTo(3.0d, TestConstants.DOUBLE_VALUE_TOLERANCE); precision = precisionMatrix.getEntry(1, 1); assertThat(precision).isEqualTo(5.5d, TestConstants.DOUBLE_VALUE_TOLERANCE); // scaled mean vector: h_x - K_xy * y RealVector scaledMeanVector = reducedVector.getScaledMeanVector(); assertThat(scaledMeanVector.getDimension()).isEqualTo(2); double meanValue = scaledMeanVector.getEntry(0); assertThat(meanValue).isEqualTo(-5.5d, TestConstants.DOUBLE_VALUE_TOLERANCE); meanValue = scaledMeanVector.getEntry(1); assertThat(meanValue).isEqualTo(-23.5d, TestConstants.DOUBLE_VALUE_TOLERANCE); // normalization constant: g + h_y * y - 0.5 * (y * K_yy * y) // 8.5 + 7.5 - 9,375 assertThat(reducedVector.getNormalizationConstant()).isEqualTo(6.625d, TestConstants.DOUBLE_VALUE_TOLERANCE); }
From source file:com.github.thorbenlindhauer.factor.CanonicalGaussianFactor.java
protected RealVector getSubVector(RealVector vector, int[] entriesToKeep) { RealVector subVector = new ArrayRealVector(entriesToKeep.length); for (int i = 0; i < entriesToKeep.length; i++) { subVector.setEntry(i, vector.getEntry(entriesToKeep[i])); }/*from w w w . j a v a2 s . c o m*/ return subVector; }
From source file:com.datumbox.framework.core.common.dataobjects.MapRealVector.java
/** {@inheritDoc} */ @Override/*from w w w. jav a 2 s . c om*/ public RealMatrix outerProduct(RealVector v) { final int m = this.getDimension(); final int n = v.getDimension(); final RealMatrix product; if (m > 1000000) { //use only in big values product = new MapRealMatrix(m, n); } else { product = new OpenMapRealMatrix(m, n); } for (int i = 0; i < m; i++) { for (int j = 0; j < n; j++) { product.setEntry(i, j, this.getEntry(i) * v.getEntry(j)); } } return product; }
From source file:edu.stanford.cfuller.imageanalysistools.clustering.ObjectClustering.java
/** * Applies basic clustering to an Image with objects. * * This will use the long-range gaussian filtering approach to assign clusters; objects sufficiently near to each other will be smeared into a single object and assigned to the same cluster. * * @param input An Image mask labeled such that each object in the Image is assigned a unique nonzero greylevel value. These should start at 1 and be consecutive. * @param original The original image (not currently used... this is here to maintain the interface with a previous version that used this image) * @param gaussianFiltered The mask with a long range Gaussian filter applied (as from {@link #gaussianFilterMask}). This is an optional parameter; * input null to have this automatically generated by the method. This parameter is * chiefly useful to save computation time when running the clutering multiple times. * This will be modified, so if planning to reuse the Gaussian filtered image, pass in a copy. *///from w w w . j a v a2 s . c om public static Image doBasicClustering(WritableImage input, Image original, Image gaussianFiltered) { RelabelFilter rlf = new RelabelFilter(); LabelFilter lf = new LabelFilter(); MaskFilter mf = new MaskFilter(); mf.setReferenceImage(input); Histogram h_individualCentromeres = new Histogram(input); WritableImage origCopy = null; if (gaussianFiltered == null) { origCopy = gaussianFilterMask(input).getWritableInstance(); } else { origCopy = gaussianFiltered.getWritableInstance(); } lf.apply(origCopy); WritableImage mapped = ImageFactory.createWritable(origCopy); Histogram h_mapped_0 = new Histogram(origCopy); //first, find the centroid of each cluster org.apache.commons.math3.linear.RealVector centroids_x = new ArrayRealVector(h_mapped_0.getMaxValue() + 1); org.apache.commons.math3.linear.RealVector centroids_y = new ArrayRealVector(h_mapped_0.getMaxValue() + 1); org.apache.commons.math3.linear.RealVector counts = new ArrayRealVector(h_mapped_0.getMaxValue() + 1); centroids_x.mapMultiplyToSelf(0.0); centroids_y.mapMultiplyToSelf(0.0); counts.mapMultiplyToSelf(0.0); for (ImageCoordinate i : origCopy) { if (origCopy.getValue(i) > 0) { int value = (int) origCopy.getValue(i); centroids_x.setEntry(value, centroids_x.getEntry(value) + i.get(ImageCoordinate.X)); centroids_y.setEntry(value, centroids_y.getEntry(value) + i.get(ImageCoordinate.Y)); counts.setEntry(value, counts.getEntry(value) + 1); } } for (int i = 0; i < counts.getDimension(); i++) { if (counts.getEntry(i) == 0) { counts.setEntry(i, 1); centroids_x.setEntry(i, -1 * origCopy.getDimensionSizes().get(ImageCoordinate.X)); centroids_y.setEntry(i, -1 * origCopy.getDimensionSizes().get(ImageCoordinate.Y)); } centroids_x.setEntry(i, centroids_x.getEntry(i) / counts.getEntry(i)); centroids_y.setEntry(i, centroids_y.getEntry(i) / counts.getEntry(i)); } for (ImageCoordinate i : origCopy) { if (mapped.getValue(i) > 0 || input.getValue(i) == 0) continue; double minDistance = Double.MAX_VALUE; int minIndex = 0; for (int j = 0; j < centroids_x.getDimension(); j++) { double dist = Math.hypot(centroids_x.getEntry(j) - i.get(ImageCoordinate.X), centroids_y.getEntry(j) - i.get(ImageCoordinate.Y)); if (dist < minDistance) { minDistance = dist; minIndex = j; } } mapped.setValue(i, minIndex); } int[] centromereAssignments = new int[h_individualCentromeres.getMaxValue() + 1]; java.util.Arrays.fill(centromereAssignments, 0); for (ImageCoordinate i : mapped) { if (input.getValue(i) > 0) { int value = (int) input.getValue(i); if (centromereAssignments[value] > 0) { mapped.setValue(i, centromereAssignments[value]); } else { centromereAssignments[value] = (int) mapped.getValue(i); } } } mf.apply(mapped); origCopy.copy(mapped); mf.setReferenceImage(origCopy); mf.apply(input); rlf.apply(input); rlf.apply(origCopy); return origCopy; }
From source file:com.datumbox.framework.core.machinelearning.regression.MatrixLinearRegression.java
/** {@inheritDoc} */ @Override//from w ww . java 2s.com protected void _fit(Dataframe trainingData) { ModelParameters modelParameters = knowledgeBase.getModelParameters(); int n = trainingData.size(); int d = trainingData.xColumnSize(); Map<Object, Double> thitas = modelParameters.getThitas(); Map<Object, Integer> featureIds = modelParameters.getFeatureIds(); //Map<Integer, Integer> recordIdsReference = null; DataframeMatrix matrixDataset = DataframeMatrix.newInstance(trainingData, true, null, featureIds); RealVector Y = matrixDataset.getY(); RealMatrix X = matrixDataset.getX(); //(X'X)^-1 RealMatrix Xt = X.transpose(); LUDecomposition lud = new LUDecomposition(Xt.multiply(X)); //W = (X'X)^-1 * X'Y RealMatrix XtXinv = lud.getSolver().getInverse(); RealVector coefficients = XtXinv.multiply(Xt).operate(Y); // instead of matrix inversion calculate (X'X) * W = X'Y //RealVector coefficients = lud.getSolver().solve(Xt.operate(Y)); //lud =null; //Xt = null; //put the features coefficients in the thita map thitas.put(Dataframe.COLUMN_NAME_CONSTANT, coefficients.getEntry(0)); for (Map.Entry<Object, Integer> entry : featureIds.entrySet()) { Object feature = entry.getKey(); Integer featureId = entry.getValue(); thitas.put(feature, coefficients.getEntry(featureId)); } //get the predictions and subtact the Y vector. Sum the squared differences to get the error double SSE = 0.0; for (double v : X.operate(coefficients).subtract(Y).toArray()) { SSE += v * v; } //Y = null; //standard error matrix double MSE = SSE / (n - (d + 1)); //mean square error = SSE / dfResidual RealMatrix SE = XtXinv.scalarMultiply(MSE); //XtXinv = null; //creating a flipped map of ids to features Map<Integer, Object> idsFeatures = PHPMethods.array_flip(featureIds); Map<Object, Double> pvalues = new HashMap<>(); //This is not small, but it does not make sense to store it in the storage for (int i = 0; i < (d + 1); ++i) { double error = SE.getEntry(i, i); Object feature = idsFeatures.get(i); if (error <= 0.0) { //double tstat = Double.MAX_VALUE; pvalues.put(feature, 0.0); } else { double tstat = coefficients.getEntry(i) / Math.sqrt(error); pvalues.put(feature, 1.0 - ContinuousDistributions.studentsCdf(tstat, n - (d + 1))); //n-d degrees of freedom } } //SE=null; //coefficients=null; //idsFeatures=null; //matrixDataset = null; modelParameters.setFeaturePvalues(pvalues); }
From source file:edu.stanford.cfuller.imageanalysistools.fitting.BisquareLinearFit.java
/** * Performs a weighted least squares fit with supplied weights to the supplied data. * * @param indVarValues A RealVector containing the values of the independent variable. * @param depVarValues A RealVector containing the values of the dependent variable. * @param weights A RealVector containing the weights for the data points. * @return a RealVector containing two elements: the slope of the fit and the y-intercept of the fit. *//* ww w . j a va2s . c om*/ public RealVector wlsFit(RealVector indVarValues, RealVector depVarValues, RealVector weights) { //initial guess for the fit: unweighted regression. SimpleRegression unweighted = new SimpleRegression(!this.noIntercept); for (int i = 0; i < indVarValues.getDimension(); i++) { unweighted.addData(indVarValues.getEntry(i), depVarValues.getEntry(i)); } RealVector parameters = null; if (this.noIntercept) { parameters = new ArrayRealVector(1, 0.0); } else { parameters = new ArrayRealVector(2, 0.0); } parameters.setEntry(0, unweighted.getSlope()); if (!this.noIntercept) { parameters.setEntry(1, unweighted.getIntercept()); } NelderMeadMinimizer nmm = new NelderMeadMinimizer(1e-12); WlsObjectiveFunction wof = new WlsObjectiveFunction(); wof.setIndVar(indVarValues); wof.setDepVar(depVarValues); wof.setWeights(weights); wof.setShouldFitIntercept(!this.noIntercept); parameters = nmm.optimize(wof, parameters); if (this.noIntercept) { RealVector output = new ArrayRealVector(2, 0.0); output.setEntry(0, parameters.getEntry(0)); return output; } return parameters; }