List of usage examples for org.apache.commons.math3.linear ArrayRealVector ArrayRealVector
public ArrayRealVector(double[] v1, double[] v2)
From source file:edu.stanford.cfuller.imageanalysistools.fitting.NelderMeadMinimizer.java
/** * Runs the minimization of the specified function starting from an initial simplex. * @param f The ObjectiveFunction to be minimized. * @param initialSimplex The initial simplex to use to start optimization, as might be returned from {@link #generateInitialSimplex} * @return The parameters at the function minimum in the same order as specified for each point on the simplex. *///from www .j av a 2 s . co m public RealVector optimize(ObjectiveFunction f, RealMatrix initialSimplex) { RealMatrix currentSimplex = initialSimplex.copy(); double currTolVal = 1.0e6; RealVector values = new ArrayRealVector(initialSimplex.getRowDimension(), 0.0); RealVector centerOfMass = new ArrayRealVector(initialSimplex.getColumnDimension(), 0.0); boolean shouldEvaluate = false; long iterCounter = 0; while (Math.abs(currTolVal) > this.tol) { int maxIndex = 0; int minIndex = 0; double maxValue = -1.0 * Double.MAX_VALUE; double minValue = Double.MAX_VALUE; double secondMaxValue = -1.0 * Double.MAX_VALUE; centerOfMass.mapMultiplyToSelf(0.0); if (shouldEvaluate) { for (int i = 0; i < currentSimplex.getRowDimension(); i++) { RealVector currRow = currentSimplex.getRowVector(i); values.setEntry(i, f.evaluate(currRow)); } } for (int i = 0; i < currentSimplex.getRowDimension(); i++) { double currValue = values.getEntry(i); if (currValue < minValue) { minValue = currValue; minIndex = i; } if (currValue > maxValue) { secondMaxValue = maxValue; maxValue = currValue; maxIndex = i; } else if (currValue > secondMaxValue) { secondMaxValue = currValue; } } for (int i = 0; i < currentSimplex.getRowDimension(); i++) { if (i == maxIndex) continue; centerOfMass = centerOfMass.add(currentSimplex.getRowVector(i)); } centerOfMass.mapDivideToSelf(currentSimplex.getRowDimension() - 1); RealVector oldPoint = currentSimplex.getRowVector(maxIndex); RealVector newPoint = centerOfMass.subtract(oldPoint).mapMultiplyToSelf(a).add(centerOfMass); // newpoint = COM + a*(COM-oldpoint) double newValue = f.evaluate(newPoint); if (newValue < secondMaxValue) { // success if (newValue < minValue) { // best found so far //expansion RealVector expPoint = centerOfMass.subtract(oldPoint).mapMultiplyToSelf(g).add(centerOfMass); double expValue = f.evaluate(expPoint); if (expValue < newValue) { currentSimplex.setRowVector(maxIndex, expPoint); currTolVal = 2.0 * (expValue - maxValue) / (1.0e-20 + expValue + maxValue); values.setEntry(maxIndex, expValue); shouldEvaluate = false; continue; } } //reflection currentSimplex.setRowVector(maxIndex, newPoint); currTolVal = 2.0 * (newValue - maxValue) / (1.0e-20 + newValue + maxValue); values.setEntry(maxIndex, newValue); shouldEvaluate = false; continue; } //contraction RealVector conPoint = centerOfMass.subtract(oldPoint).mapMultiplyToSelf(r).add(centerOfMass); double conValue = f.evaluate(conPoint); if (conValue < maxValue) { currentSimplex.setRowVector(maxIndex, conPoint); currTolVal = 2.0 * (conValue - maxValue) / (1.0e-20 + conValue + maxValue); values.setEntry(maxIndex, conValue); shouldEvaluate = false; continue; } //reduction for (int i = 0; i < currentSimplex.getRowDimension(); i++) { if (i == minIndex) continue; currentSimplex.setRowVector(i, currentSimplex.getRowVector(i).subtract(currentSimplex.getRowVector(minIndex)) .mapMultiplyToSelf(s).add(currentSimplex.getRowVector(minIndex))); } double redValue = f.evaluate(currentSimplex.getRowVector(maxIndex)); currTolVal = 2.0 * (redValue - maxValue) / (1.0e-20 + redValue + maxValue); shouldEvaluate = true; if (iterCounter++ > 100000) { System.out.println("stalled? tol: " + currTolVal + " minValue: " + minValue); } } double minValue = Double.MAX_VALUE; RealVector minVector = null; for (int i = 0; i < currentSimplex.getRowDimension(); i++) { values.setEntry(i, f.evaluate(currentSimplex.getRowVector(i))); if (values.getEntry(i) < minValue) { minValue = values.getEntry(i); minVector = currentSimplex.getRowVector(i); } } return minVector; }
From source file:edu.stanford.cfuller.colocalization3d.correction.PositionCorrector.java
/** * Creates a correction from a set of objects whose positions should be the same in each channel. * * @param imageObjects A Vector containing all the ImageObjects to be used for the correction * or in the order it appears in a multiwavelength image file. * @return A Correction object that can be used to correct the positions of other objects based upon the standards provided. *///from www .j a v a2 s.c o m public Correction getCorrection(java.util.List<ImageObject> imageObjects) { int referenceChannel = this.parameters.getIntValueForKey(REF_CH_PARAM); int channelToCorrect = this.parameters.getIntValueForKey(CORR_CH_PARAM); if (!this.parameters.hasKeyAndTrue(DET_CORR_PARAM)) { try { return Correction.readFromDisk(FileUtils.getCorrectionFilename(this.parameters)); } catch (java.io.IOException e) { java.util.logging.Logger .getLogger(edu.stanford.cfuller.colocalization3d.Colocalization3DMain.LOGGER_NAME) .severe("Exception encountered while reading correction from disk: "); e.printStackTrace(); } catch (ClassNotFoundException e) { java.util.logging.Logger .getLogger(edu.stanford.cfuller.colocalization3d.Colocalization3DMain.LOGGER_NAME) .severe("Exception encountered while reading correction from disk: "); e.printStackTrace(); } return null; } int numberOfPointsToFit = this.parameters.getIntValueForKey(NUM_POINT_PARAM); RealMatrix correctionX = new Array2DRowRealMatrix(imageObjects.size(), numberOfCorrectionParameters); RealMatrix correctionY = new Array2DRowRealMatrix(imageObjects.size(), numberOfCorrectionParameters); RealMatrix correctionZ = new Array2DRowRealMatrix(imageObjects.size(), numberOfCorrectionParameters); RealVector distanceCutoffs = new ArrayRealVector(imageObjects.size(), 0.0); RealVector ones = new ArrayRealVector(numberOfPointsToFit, 1.0); RealVector distancesToObjects = new ArrayRealVector(imageObjects.size(), 0.0); RealMatrix allCorrectionParametersMatrix = new Array2DRowRealMatrix(numberOfPointsToFit, numberOfCorrectionParameters); for (int i = 0; i < imageObjects.size(); i++) { RealVector ithPos = imageObjects.get(i).getPositionForChannel(referenceChannel); for (int j = 0; j < imageObjects.size(); j++) { double d = imageObjects.get(j).getPositionForChannel(referenceChannel).subtract(ithPos).getNorm(); distancesToObjects.setEntry(j, d); } //the sorting becomes a bottleneck once the number of points gets large //reverse comparator so we can use the priority queue and get the max element at the head Comparator<Double> cdReverse = new Comparator<Double>() { public int compare(Double o1, Double o2) { if (o1.equals(o2)) return 0; if (o1 > o2) return -1; return 1; } }; PriorityQueue<Double> pq = new PriorityQueue<Double>(numberOfPointsToFit + 2, cdReverse); double maxElement = Double.MAX_VALUE; for (int p = 0; p < numberOfPointsToFit + 1; p++) { pq.add(distancesToObjects.getEntry(p)); } maxElement = pq.peek(); for (int p = numberOfPointsToFit + 1; p < distancesToObjects.getDimension(); p++) { double value = distancesToObjects.getEntry(p); if (value < maxElement) { pq.poll(); pq.add(value); maxElement = pq.peek(); } } double firstExclude = pq.poll(); double lastDist = pq.poll(); double distanceCutoff = (lastDist + firstExclude) / 2.0; distanceCutoffs.setEntry(i, distanceCutoff); RealVector xPositionsToFit = new ArrayRealVector(numberOfPointsToFit, 0.0); RealVector yPositionsToFit = new ArrayRealVector(numberOfPointsToFit, 0.0); RealVector zPositionsToFit = new ArrayRealVector(numberOfPointsToFit, 0.0); RealMatrix differencesToFit = new Array2DRowRealMatrix(numberOfPointsToFit, imageObjects.get(0).getPositionForChannel(referenceChannel).getDimension()); int toFitCounter = 0; for (int j = 0; j < imageObjects.size(); j++) { if (distancesToObjects.getEntry(j) < distanceCutoff) { xPositionsToFit.setEntry(toFitCounter, imageObjects.get(j).getPositionForChannel(referenceChannel).getEntry(0)); yPositionsToFit.setEntry(toFitCounter, imageObjects.get(j).getPositionForChannel(referenceChannel).getEntry(1)); zPositionsToFit.setEntry(toFitCounter, imageObjects.get(j).getPositionForChannel(referenceChannel).getEntry(2)); differencesToFit.setRowVector(toFitCounter, imageObjects.get(j) .getVectorDifferenceBetweenChannels(referenceChannel, channelToCorrect)); toFitCounter++; } } RealVector x = xPositionsToFit.mapSubtractToSelf(ithPos.getEntry(0)); RealVector y = yPositionsToFit.mapSubtractToSelf(ithPos.getEntry(1)); allCorrectionParametersMatrix.setColumnVector(0, ones); allCorrectionParametersMatrix.setColumnVector(1, x); allCorrectionParametersMatrix.setColumnVector(2, y); allCorrectionParametersMatrix.setColumnVector(3, x.map(new Power(2))); allCorrectionParametersMatrix.setColumnVector(4, y.map(new Power(2))); allCorrectionParametersMatrix.setColumnVector(5, x.ebeMultiply(y)); DecompositionSolver solver = (new QRDecomposition(allCorrectionParametersMatrix)).getSolver(); RealVector cX = solver.solve(differencesToFit.getColumnVector(0)); RealVector cY = solver.solve(differencesToFit.getColumnVector(1)); RealVector cZ = solver.solve(differencesToFit.getColumnVector(2)); correctionX.setRowVector(i, cX); correctionY.setRowVector(i, cY); correctionZ.setRowVector(i, cZ); } Correction c = new Correction(correctionX, correctionY, correctionZ, distanceCutoffs, imageObjects, referenceChannel, channelToCorrect); return c; }
From source file:gamlss.distributions.GA.java
/** First derivative dlds = dl/dsigma, * where l - log-likelihood function.// w ww . ja v a 2s . c om * @param y - vector of values of response variable * @return a vector of First derivative dlds = dl/dsigma */ private ArrayRealVector dlds(final ArrayRealVector y) { //function(y,mu,sigma) (2/sigma^3)*((y/mu)-log(y) //+log(mu)+log(sigma^2)-1+digamma(1/(sigma^2))), double[] out = new double[size]; for (int i = 0; i < size; i++) { out[i] = (2 / (tempV2.getEntry(i) * sigmaV.getEntry(i))) * ((y.getEntry(i) / muV.getEntry(i)) - FastMath.log(y.getEntry(i)) + FastMath.log(muV.getEntry(i)) + FastMath.log(tempV2.getEntry(i)) - 1 + Gamma.digamma(1 / (tempV2.getEntry(i)))); } return new ArrayRealVector(out, false); }
From source file:edu.tum.cs.vis.model.util.algorithm.ACCUM.java
/** * Calculate curvature for vertices of triangle * //from www .j a v a 2s . c om * @param curvatures * vertex curvature mapping * @param tri * triangle to calculate curvature for */ static void calculateCurvatureForTriangle(HashMap<Vertex, Curvature> curvatures, Triangle tri) { // Edges Vector3f e[] = tri.getEdges(); // N-T-B coordinate system per face Vector3f t = new Vector3f(e[0]); t.normalize(); Vector3f n = new Vector3f(); n.cross(e[0], e[1]); Vector3f b = new Vector3f(); b.cross(n, t); b.normalize(); // Estimate curvature based on variation of normals // along edges float m[] = { 0, 0, 0 }; double w[][] = { { 0, 0, 0 }, { 0, 0, 0 }, { 0, 0, 0 } }; for (int j = 0; j < 3; j++) { float u = e[j].dot(t); float v = e[j].dot(b); w[0][0] += u * u; w[0][1] += u * v; // w[1][1] += v*v + u*u; // w[1][2] += u*v; w[2][2] += v * v; Vector3f dn = new Vector3f(tri.getPosition()[(j + 2) % 3].getNormalVector()); dn.sub(tri.getPosition()[(j + 1) % 3].getNormalVector()); float dnu = dn.dot(t); float dnv = dn.dot(b); m[0] += dnu * u; m[1] += dnu * v + dnv * u; m[2] += dnv * v; } w[1][1] = w[0][0] + w[2][2]; w[1][2] = w[0][1]; RealMatrix coefficients = new Array2DRowRealMatrix(w, false); DecompositionSolver solver = new LUDecomposition(coefficients).getSolver(); if (!solver.isNonSingular()) { return; } RealVector constants = new ArrayRealVector(new double[] { m[0], m[1], m[2] }, false); RealVector solution = solver.solve(constants); m[0] = (float) solution.getEntry(0); m[1] = (float) solution.getEntry(1); m[2] = (float) solution.getEntry(2); // Push it back out to the vertices for (int j = 0; j < 3; j++) { Vertex vj = tri.getPosition()[j]; float c1, c12, c2; float ret[] = proj_curv(t, b, m[0], m[1], m[2], curvatures.get(vj).getPrincipleDirectionMax(), curvatures.get(vj).getPrincipleDirectionMin()); c1 = ret[0]; c12 = ret[1]; c2 = ret[2]; Curvature c = curvatures.get(vj); double wt; if (j == 0) wt = tri.getCornerarea().x / vj.getPointarea(); else if (j == 1) wt = tri.getCornerarea().y / vj.getPointarea(); else wt = tri.getCornerarea().z / vj.getPointarea(); synchronized (c) { c.setCurvatureMax((float) (c.getCurvatureMax() + wt * c1)); c.setCurvatureMinMax((float) (c.getCurvatureMinMax() + wt * c12)); c.setCurvatureMin((float) (c.getCurvatureMin() + wt * c2)); } } }
From source file:edu.oregonstate.eecs.mcplan.ml.ConstrainedKMeans.java
private boolean mstep() { final double[][] mu_prime = new double[k][d]; final int[] nh = new int[k]; for (int i = 0; i < N; ++i) { final int h = assignments_[i]; nh[h] += 1;//from www .jav a2s . c o m final RealVector xi = X_.get(i); // mu_prime[h] = mu_prime[h].add( xi.subtract( mu_prime[h] ).mapDivideToSelf( i + 1 ) ); for (int j = 0; j < d; ++j) { // Incremental mean calculation mu_prime[h][j] += (xi.getEntry(j) - mu_prime[h][j]) / nh[h]; } } double Jprime = 0.0; for (int i = 0; i < N; ++i) { final int h = assignments_[i]; Jprime += J(i, h); } double delta = 0.0; for (int h = 0; h < k; ++h) { final RealVector v = new ArrayRealVector(mu_prime[h], false); // No copy // System.out.println( mu_[h] ); // System.out.println( v ); final RealVector diff = mu_[h].subtract(v); delta += diff.getNorm(); // delta += mu_[h].getDistance( v ); // delta += distance( mu_[h], v ); mu_[h] = v; } // System.out.println( "\tdelta = " + delta ); final double deltaJ = J_ - Jprime; // System.out.println( "\tdeltaJ = " + deltaJ ); J_ = Jprime; return deltaJ > convergence_tolerance_; }
From source file:gamlss.distributions.TF2.java
/** Set sigma1 array. * @param y - response variable//from w w w . j a v a 2 s. c o m */ private void setInterimArrays(final ArrayRealVector y) { muV = distributionParameters.get(DistributionSettings.MU); sigmaV = distributionParameters.get(DistributionSettings.SIGMA); nuV = distributionParameters.get(DistributionSettings.NU); size = y.getDimension(); double[] temp = new double[size]; for (int i = 0; i < size; i++) { //sigma1 <- (sqrt((nu-2)/nu))*sigma temp[i] = (FastMath.sqrt((nuV.getEntry(i) - 2.0) / nuV.getEntry(i))) * sigmaV.getEntry(i); } sigma1 = new ArrayRealVector(temp, false); }
From source file:gamlss.distributions.GA.java
/** Second derivative d2ldm2= (d^2l)/(dmu^2), * where l - log-likelihood function./*from ww w . ja v a2s .co m*/ * @return a vector of Second derivative d2ldm2= (d^2l)/(dmu^2) */ private ArrayRealVector d2ldm2() { //d2ldm2 = function(mu,sigma) -1/((sigma^2)*(mu^2)) double[] out = new double[size]; for (int i = 0; i < size; i++) { out[i] = -1 / (tempV2.getEntry(i) * (muV.getEntry(i) * muV.getEntry(i))); } return new ArrayRealVector(out, false); }
From source file:gamlss.distributions.NO.java
/** First derivative dldm = dl/dmu, * where l - log-likelihood function./*from w ww . ja v a 2s .com*/ * @param y - vector of values of response variable * @return a vector of first derivative dldm = dl/dmu */ private ArrayRealVector dldm(final ArrayRealVector y) { double[] out = new double[size]; for (int i = 0; i < size; i++) { //(1/sigma^2)*(y-mu) out[i] = (1 / tempV2.getEntry(i)) * (y.getEntry(i) - muV.getEntry(i)); } return new ArrayRealVector(out, false); }
From source file:gamlss.distributions.TF.java
/** First derivative dldm = dl/dmu, where l - log-likelihood function. * @param y - vector of values of response variable * @return a vector of first derivative dldm = dl/dmu *//*from w w w . j a va2 s. co m*/ public final ArrayRealVector dldm(final ArrayRealVector y) { double[] dldm = new double[size]; for (int i = 0; i < size; i++) { //dldm <- (omega*(y-mu))/s2 dldm[i] = (omega[i] * ym[i]) / s2[i]; } dsq = null; omega = null; return new ArrayRealVector(dldm, false); }
From source file:gamlss.algorithm.GlimFitCG.java
private ArrayRealVector setWV(ArrayRealVector in1, ArrayRealVector in2) { double[] out = new double[in1.getDimension()]; for (int i = 0; i < out.length; i++) { out[i] = in1.getEntry(i) + in2.getEntry(i); }//from w w w . j a va 2 s .c o m return new ArrayRealVector(out, false); }