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.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. *//*www. ja v a 2s. 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:edu.stanford.cfuller.imageanalysistools.clustering.GaussianLikelihoodObjectiveFunction.java
/** * Evaluates the function with the specified parameters. * * The parameters describe a set of gaussian generators (which are the Clusters). * * @param parameters A RealVector containing the values of all the parameters of each Gaussian, ordered so that all the parameters of a single gaussian are together, then the next gaussian, etc. * @return The negative log likelihood of having observed the ClusterObjects at their locations, given the parameters describing the Gaussian clusters. */// ww w . j ava 2s . c o m public double evaluate(RealVector parameters) { int nClusters = parameters.getDimension() / nParametersEach; //java.util.logging.Logger.getLogger("edu.stanford.cfuller.imageanalysistools").info("nClusters: " + nClusters + " abdMatrices_size: " + abdMatrices.size() + " det_dim: " + det.getDimension()); if (det.getDimension() != nClusters) { clusterProbs = new Array2DRowRealMatrix(this.objects.size(), nClusters); det = new ArrayRealVector(nClusters); pk = new ArrayRealVector(nClusters); if (abdMatrices.size() < nClusters) { int originalSize = abdMatrices.size(); for (int i = 0; i < nClusters - originalSize; i++) { abdMatrices.add(new Array2DRowRealMatrix(numDim, numDim)); } } else { abdMatrices.setSize(nClusters); } } pk.mapMultiplyToSelf(0.0); //java.util.logging.Logger.getLogger("edu.stanford.cfuller.imageanalysistools").info("nClusters: " + nClusters + " abdMatrices_size: " + abdMatrices.size() + " det_dim: " + det.getDimension()); for (int i = 0; i < nClusters; i++) { /* double ct = Math.cos(parameters.getEntry(nParametersEach*i+3)); double st = Math.sin(parameters.getEntry(nParametersEach*i+3)); double sin2t = Math.sin(2*parameters.getEntry(nParametersEach*i+3)); double a = (ct*ct/(2*parameters.getEntry(nParametersEach*i+2)) + st*st/(2*parameters.getEntry(nParametersEach*i+4))); double b = (sin2t/(4*parameters.getEntry(nParametersEach*i+4)) - sin2t/(4*parameters.getEntry(nParametersEach*i+2))); double d = (st*st/(2*parameters.getEntry(nParametersEach*i+2)) + ct*ct/(2*parameters.getEntry(nParametersEach*i+4))); */ double a = parameters.getEntry(nParametersEach * i + 2); double d = parameters.getEntry(nParametersEach * i + 4); double b = Math.sqrt(a * d) * parameters.getEntry(nParametersEach * i + 3); abdMatrices.get(i).setEntry(0, 0, a); abdMatrices.get(i).setEntry(1, 0, b); abdMatrices.get(i).setEntry(0, 1, b); abdMatrices.get(i).setEntry(1, 1, d); LUDecomposition abdLU = (new LUDecomposition(abdMatrices.get(i))); det.setEntry(i, (abdLU).getDeterminant()); //det.setEntry(i, a*d-b*b); try { abdMatrices.set(i, abdLU.getSolver().getInverse()); } catch (org.apache.commons.math3.linear.SingularMatrixException e) { return Double.MAX_VALUE; } } for (int n = 0; n < this.objects.size(); n++) { ClusterObject c = this.objects.get(n); double max = -1.0 * Double.MAX_VALUE; int maxIndex = 0; for (int k = 0; k < nClusters; k++) { mean.setEntry(0, c.getCentroid().getX() - parameters.getEntry(nParametersEach * k)); mean.setEntry(1, c.getCentroid().getY() - parameters.getEntry(nParametersEach * k + 1)); x = abdMatrices.get(k).operate(mean); double dot = x.dotProduct(mean); // java.util.logging.Logger.getLogger("edu.stanford.cfuller.imageanalysistools").info("k, n: " + k + ", " + this.objects.size()); // java.util.logging.Logger.getLogger("edu.stanford.cfuller.imageanalysistools").info("parameters: " + parameters.toString()); // java.util.logging.Logger.getLogger("edu.stanford.cfuller.imageanalysistools").info("abd matrix: " + abdMatrices.get(k).toString()); // java.util.logging.Logger.getLogger("edu.stanford.cfuller.imageanalysistools").info("det: " + det.getEntry(k)); // java.util.logging.Logger.getLogger("edu.stanford.cfuller.imageanalysistools").info("mean: " + mean.toString()); // java.util.logging.Logger.getLogger("edu.stanford.cfuller.imageanalysistools").info("dot: " + dot); double logN = negLog2PI - 0.5 * Math.log(det.getEntry(k)) - 0.5 * dot; // java.util.logging.Logger.getLogger("edu.stanford.cfuller.imageanalysistools").info("logN: " + logN); clusterProbs.setEntry(n, k, logN); if (logN > max) { max = logN; maxIndex = k; } if (Double.isInfinite(logN) || Double.isNaN(logN)) { return Double.MAX_VALUE; } } c.setMostProbableCluster(maxIndex); } for (int k = 0; k < nClusters; k++) { double tempMax = -1.0 * Double.MAX_VALUE; for (int n = 0; n < this.objects.size(); n++) { if (clusterProbs.getEntry(n, k) > tempMax) tempMax = clusterProbs.getEntry(n, k); } pk.setEntry(k, tempMax + Math.log( clusterProbs.getColumnVector(k).mapSubtract(tempMax).mapToSelf(new Exp()).getL1Norm()) - Math.log(this.objects.size())); } double pkMax = -1.0 * Double.MAX_VALUE; for (int k = 0; k < nClusters; k++) { if (pk.getEntry(k) > pkMax) pkMax = pk.getEntry(k); } double logSumPk = pkMax + Math.log(pk.mapSubtract(pkMax).mapToSelf(new Exp()).getL1Norm()); pk.mapSubtractToSelf(logSumPk); //java.util.logging.Logger.getLogger("edu.stanford.cfuller.imageanalysistools").info("pk: " + pk.toString()); double L = 0; for (int n = 0; n < this.objects.size(); n++) { RealVector toSum = clusterProbs.getRowVector(n).add(pk); double tempMax = -1.0 * Double.MAX_VALUE; for (int k = 0; k < nClusters; k++) { if (toSum.getEntry(k) > tempMax) tempMax = toSum.getEntry(k); } double pn = tempMax + Math.log(toSum.mapSubtract(tempMax).mapToSelf(new Exp()).getL1Norm()); //java.util.logging.Logger.getLogger("edu.stanford.cfuller.imageanalysistools").info("pn: " + pn); L += pn; } return -1.0 * L; }
From source file:edu.oregonstate.eecs.mcplan.abstraction.Experiments.java
private static void writeClustering(final MetricConstrainedKMeans kmeans, final File root, final int iter) throws FileNotFoundException { Csv.write(new PrintStream(new File(root, "M" + iter + ".csv")), kmeans.metric); {/*from www . j a v a2s. co m*/ final Csv.Writer writer = new Csv.Writer(new PrintStream(new File(root, "mu" + iter + ".csv"))); for (int i = 0; i < kmeans.d; ++i) { for (int j = 0; j < kmeans.k; ++j) { writer.cell(kmeans.mu()[j].getEntry(i)); } writer.newline(); } } // Lt.operate( x ) maps x to the space defined by the metric final RealMatrix Lt = new CholeskyDecomposition(kmeans.metric).getLT(); { final Csv.Writer writer = new Csv.Writer(new PrintStream(new File(root, "X" + iter + ".csv"))); writer.cell("cluster").cell("label"); for (int i = 0; i < kmeans.metric.getColumnDimension(); ++i) { writer.cell("x" + i); } for (int i = 0; i < kmeans.metric.getColumnDimension(); ++i) { writer.cell("Ax" + i); } writer.newline(); for (int cluster = 0; cluster < kmeans.k; ++cluster) { for (int i = 0; i < kmeans.N; ++i) { if (kmeans.assignments()[i] == cluster) { writer.cell(cluster); final RealVector phi = kmeans.X_.get(i); //Phi.get( i ); writer.cell("?"); // TODO: write label for (int j = 0; j < phi.getDimension(); ++j) { writer.cell(phi.getEntry(j)); } final RealVector trans = Lt.operate(phi); for (int j = 0; j < trans.getDimension(); ++j) { writer.cell(trans.getEntry(j)); } writer.newline(); } } } } }
From source file:edu.stanford.cfuller.imageanalysistools.fitting.DifferentialEvolutionMinimizer.java
/** * Performs a minimization of a function starting with a given population. * /*from w w w . j a va 2s.c o m*/ * @param population The population of parameters to start from, one population entry per row, one parameter per column. * @param f The function to be minimized. * @param parameterLowerBounds The lower bounds of each parameter. This must have the same size as the column dimension of the population. Generated parameter values less than these values will be discarded. * @param parameterUpperBounds The upper bounds of each paraemter. This must have the same size as the column dimension of the population. Generated parameter values greater than these values will be discarded. * @param populationSize The size of the population of parameters sets. This must be equal to the row dimension of the population. * @param scaleFactor Factor controlling the scale of crossed over entries during crossover events. * @param maxIterations The maximum number of iterations to allow before returning a result. * @param crossoverFrequency The frequency of crossover from 0 to 1. At any given parameter, this is the probability of initiating a crossover as well as the probability of ending one after it has started. * @param tol Relative function value tolerance controlling termination; if the maximal and minimal population values differ by less than this factor times the maximal value, optimization will terminate. * @return The parameter values at the minimal function value found. */ public RealVector minimizeWithPopulation(RealMatrix population, ObjectiveFunction f, RealVector parameterLowerBounds, RealVector parameterUpperBounds, int populationSize, double scaleFactor, int maxIterations, double crossoverFrequency, double tol) { int numberOfParameters = parameterLowerBounds.getDimension(); double currMinValue = Double.MAX_VALUE; double currMaxValue = -1.0 * Double.MAX_VALUE; int iterationCounter = maxIterations; double mutationProb = 0.01; // int totalIterations =0; RealVector values = new ArrayRealVector(populationSize); boolean[] valuesChanged = new boolean[populationSize]; java.util.Arrays.fill(valuesChanged, true); computeValues(f, population, values, valuesChanged); RealVector newVec = new ArrayRealVector(numberOfParameters); RealMatrix newPopulation = new Array2DRowRealMatrix(populationSize, numberOfParameters); while (iterationCounter > 0) { for (int i = 0; i < populationSize; i++) { int i1 = RandomGenerator.getGenerator().randInt(populationSize); int i2 = RandomGenerator.getGenerator().randInt(populationSize); int i3 = RandomGenerator.getGenerator().randInt(populationSize); newVec.mapMultiplyToSelf(0.0); boolean inBounds = true; boolean isCrossingOver = false; for (int j = 0; j < numberOfParameters; j++) { if ((RandomGenerator.rand() < crossoverFrequency) ^ isCrossingOver) { if (!isCrossingOver) { isCrossingOver = true; } newVec.setEntry(j, scaleFactor * (population.getEntry(i2, j) - population.getEntry(i1, j)) + population.getEntry(i3, j)); } else { if (isCrossingOver) { isCrossingOver = false; } newVec.setEntry(j, population.getEntry(i, j)); } //random 10% range +/- mutation if ((RandomGenerator.rand() < mutationProb)) { double magnitude = 0.2 * (parameterUpperBounds.getEntry(j) - parameterLowerBounds.getEntry(j)); newVec.setEntry(j, newVec.getEntry(j) + (RandomGenerator.rand() - 0.5) * magnitude); } if (newVec.getEntry(j) < parameterLowerBounds.getEntry(j) || newVec.getEntry(j) > parameterUpperBounds.getEntry(j)) { inBounds = false; } } double functionValue = Double.MAX_VALUE; if (inBounds) functionValue = f.evaluate(newVec); //if (inBounds) System.out.printf("in bounds candidate value: %1.2f old value: %1.2f \n", functionValue, values.getEntry(i)); if (functionValue < values.getEntry(i)) { newPopulation.setRowVector(i, newVec); valuesChanged[i] = true; values.setEntry(i, functionValue); } else { newPopulation.setRowVector(i, population.getRowVector(i)); valuesChanged[i] = false; } } population = newPopulation; double tempMinValue = Double.MAX_VALUE; double tempMaxValue = -1.0 * Double.MAX_VALUE; // double averageValue = 0; for (int i = 0; i < values.getDimension(); i++) { double value = values.getEntry(i); // averageValue += value; if (value < tempMinValue) { tempMinValue = value; } if (value > tempMaxValue) { tempMaxValue = value; } } // averageValue /= values.getDimension(); currMinValue = tempMinValue; currMaxValue = tempMaxValue; //LoggingUtilities.getLogger().info("Iteration counter: " + Integer.toString(totalIterations) + " best score: " + currMinValue + " worst score: " + currMaxValue + " average score: " + averageValue); if (Math.abs(currMaxValue - currMinValue) < Math.abs(tol * currMaxValue) + Math.abs(tol * currMinValue)) { iterationCounter--; } else { iterationCounter = 1; } // totalIterations++; } for (int i = 0; i < populationSize; i++) { valuesChanged[i] = true; } computeValues(f, population, values, valuesChanged); double tempMinValue = Double.MAX_VALUE; int tempMinIndex = 0; for (int i = 0; i < populationSize; i++) { if (values.getEntry(i) < tempMinValue) { tempMinValue = values.getEntry(i); tempMinIndex = i; } } RealVector output = new ArrayRealVector(numberOfParameters); for (int i = 0; i < numberOfParameters; i++) { output.setEntry(i, population.getEntry(tempMinIndex, i)); } return output; }
From source file:monitor.processing.OLD.Drawing3D.java
public void drawTree(BaseNode n, RealVector v) { Color c = Color.WHITE; if (n instanceof TransformNode) { //v = ((TransformNode) n).getTrasformMatrix().getColumnVector(3); v = ((TransformNode) n).getTrasformMatrix().operate(v); //stroke(0, 0, 255); // line((float) bef[0] * scale, (float) bef[1] * scale, (float) bef[2] * scale, (float) att[0] * scale, (float) att[1] * scale, (float) att[2] * scale); c = Color.BLACK;//from www . ja v a 2 s . c o m } else if (n instanceof StaticMesh) { StaticMesh sm = (StaticMesh) n; if (sm.isTransparent()) { c = Color.BLUE; } else if (sm.isVisible()) { // System.out.println(n.getInfo()); // System.out.println(v); c = Color.RED; if (sm.getModel().contains("naohead") || sm.getModel().contains("naobody")) { c = Color.PINK; } } else { c = Color.CYAN; } } else if (n instanceof StaticMeshNode) { StaticMeshNode smn = (StaticMeshNode) n; if (smn.isTransparent()) { c = Color.GREEN; // cameraSceneX = (float)v.getEntry(0); // cameraSceneY = (float)v.getEntry(1); // cameraSceneZ = (float)v.getEntry(2); } else if (smn.isVisible()) { c = Color.YELLOW; } else { c = Color.MAGENTA; } } // if (n.getInfo().contains("naohead")){ // c = Color.CYAN; // pushMatrix(); // fill(c.getRed(), c.getGreen(), c.getBlue()); // // noStroke(); // translate((float) (scale * v.getEntry(0)), (float) (-scale * v.getEntry(1)), (float) (scale * v.getEntry(2))); // box(50); // popMatrix(); // } pushMatrix(); fill(c.getRed(), c.getGreen(), c.getBlue()); noStroke(); translate((float) (scale * v.getEntry(0)), (float) (-scale * v.getEntry(1)), (float) (scale * v.getEntry(2))); box(0.05f * scale); popMatrix(); if (n.getAddress() == selected) { path.add(v); pushMatrix(); fill(0, 255, 255, selectedAlpha); selectedAlpha += (selectedAlpha >= 255) ? -255 : 10; noStroke(); translate((float) (scale * v.getEntry(0)), (float) (-scale * v.getEntry(1)), (float) (scale * v.getEntry(2))); pushMatrix(); //rotateX(t); rotateZ(selectedTheta); selectedTheta += 1 / 8f; box(0.1f * scale); popMatrix(); // info rotateZ(PI / 2); rotateX(-PI / 2); translate(0.1f * scale, -0.1f * scale); fill(0); textSize(0.1f * scale); text(n.getInfo(), 0, 0, 0); popMatrix(); RealVector ant = null; int o = 0; for (RealVector w : path) { if (ant != null) { stroke(o, 255 - o, 0, 200); // System.out.println(ant); // System.out.println(w); line((float) ant.getEntry(0) * scale, (float) ant.getEntry(1) * -scale, (float) ant.getEntry(2) * scale, (float) w.getEntry(0) * scale, (float) w.getEntry(1) * -scale, (float) w.getEntry(2)); o += (o >= 255) ? -200 : 10; } ant = w; pushMatrix(); stroke(255, 255, 0, 200); translate((float) (scale * w.getEntry(0)), (float) (-scale * w.getEntry(1)), (float) (scale * w.getEntry(2))); box(2); popMatrix(); } } for (int a = n.getChildren().size() - 1; a >= 0; a--) { drawTree(((ArrayList<BaseNode>) n.getChildren()).get(a), v); } }
From source file:com.winterwell.maths.stats.algorithms.KalmanFilterTest.java
@Test public void testConstantAcceleration() { // simulates a vehicle, accelerating at a constant rate (0.1 m/s) // discrete time interval double dt = 0.1d; // position measurement noise (meter) double measurementNoise = 10d; // acceleration noise (meter/sec^2) double accelNoise = 0.2d; // A = [ 1 dt ] // [ 0 1 ] RealMatrix A = new Array2DRowRealMatrix(new double[][] { { 1, dt }, { 0, 1 } }); // B = [ dt^2/2 ] // [ dt ] RealMatrix Bnull = new Array2DRowRealMatrix(new double[][] { { FastMath.pow(dt, 2d) / 2d }, { dt } }); // H = [ 1 0 ] RealMatrix H = new Array2DRowRealMatrix(new double[][] { { 1d, 0d } }); // x = [ 0 0 ] RealVector x = new ArrayRealVector(new double[] { 0, 0 }); RealMatrix tmp = new Array2DRowRealMatrix( new double[][] { { FastMath.pow(dt, 4d) / 4d, FastMath.pow(dt, 3d) / 2d }, { FastMath.pow(dt, 3d) / 2d, FastMath.pow(dt, 2d) } }); // Q = [ dt^4/4 dt^3/2 ] // [ dt^3/2 dt^2 ] RealMatrix Q = tmp.scalarMultiply(FastMath.pow(accelNoise, 2)); // P0 = [ 1 1 ] // [ 1 1 ] RealMatrix P0 = new Array2DRowRealMatrix(new double[][] { { 1, 1 }, { 1, 1 } }); // R = [ measurementNoise^2 ] RealMatrix R = new Array2DRowRealMatrix(new double[] { FastMath.pow(measurementNoise, 2) }); // constant control input, increase velocity by 0.1 m/s per cycle double uv = 0.1d * dt; RealVector u = new ArrayRealVector(new double[] { uv * uv / 2, uv }); ProcessModel pm = new DefaultProcessModel(A, Bnull, Q, x, P0); MeasurementModel mm = new DefaultMeasurementModel(H, R); KalmanFilter filter = new KalmanFilter(pm, mm); Assert.assertEquals(1, filter.getMeasurementDimension()); Assert.assertEquals(2, filter.getStateDimension()); Gaussian state = new Gaussian(mtj(x), mtj(P0)); MatrixUtils.equals(mtj(P0), state.getCovar()); // check the initial state double[] expectedInitialState = new double[] { 0.0, 0.0 }; // assertVectorEquals(expectedInitialState, filter.getStateEstimation()); RandomGenerator rand = new JDKRandomGenerator(); RealVector tmpPNoise = new ArrayRealVector(new double[] { FastMath.pow(dt, 2d) / 2d, dt }); // iterate 60 steps for (int i = 0; i < 60; i++) { state = filter.predict(state, mtj(u)); // Simulate the process RealVector pNoise = tmpPNoise.mapMultiply(accelNoise * rand.nextGaussian()); // x = A * x + B * u + pNoise x = A.operate(x).add(u).add(pNoise); // Simulate the measurement double mNoise = measurementNoise * rand.nextGaussian(); // z = H * x + m_noise RealVector z = H.operate(x).mapAdd(mNoise); state = filter.correct(state, mtj(z)); // state estimate shouldn't be larger than the measurement noise double diff = FastMath.abs(x.getEntry(0) - state.getMean().get(0)); Assert.assertTrue(Precision.compareTo(diff, measurementNoise, 1e-6) < 0); }/*from w w w. j a v a2s .c o m*/ // error covariance of the velocity should be already very low (< 0.1) Assert.assertTrue(Precision.compareTo(state.getCovar().get(1, 1), 0.1d, 1e-6) < 0); }
From source file:com.joptimizer.optimizers.LPPrimalDualMethodNetlibTest.java
private void checkSolution(LPNetlibProblem problem, LPOptimizationRequest or, LPOptimizationResponse response) throws Exception { double expectedvalue = problem.optimalValue; log.debug("expectedvalue : " + expectedvalue); double[] sol = response.getSolution(); RealVector cVector = new ArrayRealVector(or.getC().toArray()); RealVector solVector = new ArrayRealVector(sol); double value = cVector.dotProduct(solVector); log.debug("sol : " + ArrayUtils.toString(sol)); log.debug("value : " + value); //check constraints assertEquals(or.getLb().size(), sol.length); assertEquals(or.getUb().size(), sol.length); RealVector x = MatrixUtils.createRealVector(sol); //x >= lb/* w w w. jav a 2s . c o m*/ double maxLbmx = -Double.MAX_VALUE; for (int i = 0; i < or.getLb().size(); i++) { double di = Double.isNaN(or.getLb().getQuick(i)) ? -Double.MAX_VALUE : or.getLb().getQuick(i); maxLbmx = Math.max(maxLbmx, di - x.getEntry(i)); assertTrue(di <= x.getEntry(i) + or.getTolerance()); } log.debug("max(lb - x): " + maxLbmx); //x <= ub double maxXmub = -Double.MAX_VALUE; for (int i = 0; i < or.getUb().size(); i++) { double di = Double.isNaN(or.getUb().getQuick(i)) ? Double.MAX_VALUE : or.getUb().getQuick(i); maxXmub = Math.max(maxXmub, x.getEntry(i) - di); assertTrue(di + or.getTolerance() >= x.getEntry(i)); } log.debug("max(x - ub): " + maxXmub); //G.x <h if (or.getG() != null && or.getG().rows() > 0) { RealMatrix GMatrix = MatrixUtils.createRealMatrix(or.getG().toArray()); RealVector hvector = MatrixUtils.createRealVector(or.getH().toArray()); RealVector Gxh = GMatrix.operate(x).subtract(hvector); double maxGxh = -Double.MAX_VALUE; for (int i = 0; i < Gxh.getDimension(); i++) { maxGxh = Math.max(maxGxh, Gxh.getEntry(i)); assertTrue(Gxh.getEntry(i) - or.getTolerance() <= 0); } log.debug("max(G.x - h): " + maxGxh); } //A.x = b if (or.getA() != null && or.getA().rows() > 0) { RealMatrix AMatrix = MatrixUtils.createRealMatrix(or.getA().toArray()); RealVector bvector = MatrixUtils.createRealVector(or.getB().toArray()); RealVector Axb = AMatrix.operate(x).subtract(bvector); double norm = Axb.getNorm(); log.debug("||A.x -b||: " + norm); assertEquals(0., norm, or.getToleranceFeas()); } double percDelta = Math.abs((expectedvalue - value) / expectedvalue); log.debug("percDelta: " + percDelta); //assertEquals(0., percDelta, or.getTolerance()); //assertEquals(expectedvalue, value, or.getTolerance()); assertTrue(value < expectedvalue + or.getTolerance());//can even beat other optimizers! the rebel yell... }
From source file:com.itemanalysis.psychometrics.optimization.BOBYQAOptimizer.java
/** * A version of the truncated conjugate gradientAt is applied. If a line * search is restricted by a constraint, then the procedure is restarted, * the values of the variables that are at their bounds being fixed. If * the trust region boundary is reached, then further changes may be made * to D, each one being in the two dimensional space that is spanned * by the current D and the gradientAt of Q at XOPT+D, staying on the trust * region boundary. Termination occurs when the reduction in Q seems to * be close to the greatest reduction that can be achieved. * The arguments N, NPT, XPT, XOPT, GOPT, HQ, PQ, SL and SU have the same * meanings as the corresponding arguments of BOBYQB. * DELTA is the trust region radius for the present calculation, which * seeks a small value of the quadratic model within distance DELTA of * XOPT subject to the bounds on the variables. * XNEW will be set to a new vector of variables that is approximately * the one that minimizes the quadratic model within the trust region * subject to the SL and SU constraints on the variables. It satisfies * as equations the bounds that become active during the calculation. * D is the calculated trial step from XOPT, generated iteratively from an * initial value of zero. Thus XNEW is XOPT+D after the final iteration. * GNEW holds the gradientAt of the quadratic model at XOPT+D. It is updated * when D is updated./*from www.j a v a 2s . co m*/ * xbdi.get( is a working space vector. For I=1,2,...,N, the element xbdi.get((I) is * set to -1.0, 0.0, or 1.0, the value being nonzero if and only if the * I-th variable has become fixed at a bound, the bound being SL(I) or * SU(I) in the case xbdi.get((I)=-1.0 or xbdi.get((I)=1.0, respectively. This * information is accumulated during the construction of XNEW. * The arrays S, HS and HRED are also used for working space. They hold the * current search direction, and the changes in the gradientAt of Q along S * and the reduced D, respectively, where the reduced D is the same as D, * except that the components of the fixed variables are zero. * DSQ will be set to the square of the length of XNEW-XOPT. * CRVMIN is set to zero if D reaches the trust region boundary. Otherwise * it is set to the least curvature of H that occurs in the conjugate * gradientAt searches that are not restricted by any constraints. The * value CRVMIN=-1.0D0 is set, however, if all of these searches are * constrained. * @param delta * @param gnew * @param xbdi * @param s * @param hs * @param hred */ private double[] trsbox(double delta, ArrayRealVector gnew, ArrayRealVector xbdi, ArrayRealVector s, ArrayRealVector hs, ArrayRealVector hred) { printMethod(); // XXX final int n = currentBest.getDimension(); final int npt = numberOfInterpolationPoints; double dsq = Double.NaN; double crvmin = Double.NaN; // Local variables double ds; int iu; double dhd, dhs, cth, shs, sth, ssq, beta = 0, sdec, blen; int iact = -1; int nact = 0; double angt = 0, qred; int isav; double temp = 0, xsav = 0, xsum = 0, angbd = 0, dredg = 0, sredg = 0; int iterc; double resid = 0, delsq = 0, ggsav = 0, tempa = 0, tempb = 0, redmax = 0, dredsq = 0, redsav = 0, gredsq = 0, rednew = 0; int itcsav = 0; double rdprev = 0, rdnext = 0, stplen = 0, stepsq = 0; int itermax = 0; // Set some constants. // Function Body // The sign of GOPT(I) gives the sign of the change to the I-th variable // that will reduce Q from its value at XOPT. Thus xbdi.get((I) shows whether // or not to fix the I-th variable at one of its bounds initially, with // NACT being set to the number of fixed variables. D and GNEW are also // set for the first iteration. DELSQ is the upper bound on the sum of // squares of the free variables. QRED is the reduction in Q so far. iterc = 0; nact = 0; for (int i = 0; i < n; i++) { xbdi.setEntry(i, ZERO); if (trustRegionCenterOffset.getEntry(i) <= lowerDifference.getEntry(i)) { if (gradientAtTrustRegionCenter.getEntry(i) >= ZERO) { xbdi.setEntry(i, MINUS_ONE); } } else if (trustRegionCenterOffset.getEntry(i) >= upperDifference.getEntry(i)) { if (gradientAtTrustRegionCenter.getEntry(i) <= ZERO) { xbdi.setEntry(i, ONE); } } if (xbdi.getEntry(i) != ZERO) { ++nact; } trialStepPoint.setEntry(i, ZERO); gnew.setEntry(i, gradientAtTrustRegionCenter.getEntry(i)); } delsq = delta * delta; qred = ZERO; crvmin = MINUS_ONE; // Set the next search direction of the conjugate gradientAt method. It is // the steepest descent direction initially and when the iterations are // restarted because a variable has just been fixed by a bound, and of // course the components of the fixed variables are zero. ITERMAX is an // upper bound on the indices of the conjugate gradientAt iterations. int state = 20; for (;;) { switch (state) { case 20: { printState(20); // XXX beta = ZERO; } case 30: { printState(30); // XXX stepsq = ZERO; for (int i = 0; i < n; i++) { if (xbdi.getEntry(i) != ZERO) { s.setEntry(i, ZERO); } else if (beta == ZERO) { s.setEntry(i, -gnew.getEntry(i)); } else { s.setEntry(i, beta * s.getEntry(i) - gnew.getEntry(i)); } // Computing 2nd power final double d1 = s.getEntry(i); stepsq += d1 * d1; } if (stepsq == ZERO) { state = 190; break; } if (beta == ZERO) { gredsq = stepsq; itermax = iterc + n - nact; } if (gredsq * delsq <= qred * 1e-4 * qred) { state = 190; break; } // Multiply the search direction by the second derivative matrix of Q and // calculate some scalars for the choice of steplength. Then set BLEN to // the length of the the step to the trust region boundary and STPLEN to // the steplength, ignoring the simple bounds. state = 210; break; } case 50: { printState(50); // XXX resid = delsq; ds = ZERO; shs = ZERO; for (int i = 0; i < n; i++) { if (xbdi.getEntry(i) == ZERO) { // Computing 2nd power final double d1 = trialStepPoint.getEntry(i); resid -= d1 * d1; ds += s.getEntry(i) * trialStepPoint.getEntry(i); shs += s.getEntry(i) * hs.getEntry(i); } } if (resid <= ZERO) { state = 90; break; } temp = Math.sqrt(stepsq * resid + ds * ds); if (ds < ZERO) { blen = (temp - ds) / stepsq; } else { blen = resid / (temp + ds); } stplen = blen; if (shs > ZERO) { // Computing MIN stplen = Math.min(blen, gredsq / shs); } // Reduce STPLEN if necessary in order to preserve the simple bounds, // letting IACT be the index of the new constrained variable. iact = -1; for (int i = 0; i < n; i++) { if (s.getEntry(i) != ZERO) { xsum = trustRegionCenterOffset.getEntry(i) + trialStepPoint.getEntry(i); if (s.getEntry(i) > ZERO) { temp = (upperDifference.getEntry(i) - xsum) / s.getEntry(i); } else { temp = (lowerDifference.getEntry(i) - xsum) / s.getEntry(i); } if (temp < stplen) { stplen = temp; iact = i; } } } // Update CRVMIN, GNEW and D. Set SDEC to the decrease that occurs in Q. sdec = ZERO; if (stplen > ZERO) { ++iterc; temp = shs / stepsq; if (iact == -1 && temp > ZERO) { crvmin = Math.min(crvmin, temp); if (crvmin == MINUS_ONE) { crvmin = temp; } } ggsav = gredsq; gredsq = ZERO; for (int i = 0; i < n; i++) { gnew.setEntry(i, gnew.getEntry(i) + stplen * hs.getEntry(i)); if (xbdi.getEntry(i) == ZERO) { // Computing 2nd power final double d1 = gnew.getEntry(i); gredsq += d1 * d1; } trialStepPoint.setEntry(i, trialStepPoint.getEntry(i) + stplen * s.getEntry(i)); } // Computing MAX final double d1 = stplen * (ggsav - HALF * stplen * shs); sdec = Math.max(d1, ZERO); qred += sdec; } // Restart the conjugate gradientAt method if it has hit a new bound. if (iact >= 0) { ++nact; xbdi.setEntry(iact, ONE); if (s.getEntry(iact) < ZERO) { xbdi.setEntry(iact, MINUS_ONE); } // Computing 2nd power final double d1 = trialStepPoint.getEntry(iact); delsq -= d1 * d1; if (delsq <= ZERO) { state = 190; break; } state = 20; break; } // If STPLEN is less than BLEN, then either apply another conjugate // gradientAt iteration or RETURN. if (stplen < blen) { if (iterc == itermax) { state = 190; break; } if (sdec <= qred * .01) { state = 190; break; } beta = gredsq / ggsav; state = 30; break; } } case 90: { printState(90); // XXX crvmin = ZERO; // Prepare for the alternative iteration by calculating some scalars // and by multiplying the reduced D by the second derivative matrix of // Q, where S holds the reduced D in the call of GGMULT. } case 100: { printState(100); // XXX if (nact >= n - 1) { state = 190; break; } dredsq = ZERO; dredg = ZERO; gredsq = ZERO; for (int i = 0; i < n; i++) { if (xbdi.getEntry(i) == ZERO) { // Computing 2nd power double d1 = trialStepPoint.getEntry(i); dredsq += d1 * d1; dredg += trialStepPoint.getEntry(i) * gnew.getEntry(i); // Computing 2nd power d1 = gnew.getEntry(i); gredsq += d1 * d1; s.setEntry(i, trialStepPoint.getEntry(i)); } else { s.setEntry(i, ZERO); } } itcsav = iterc; state = 210; break; // Let the search direction S be a linear combination of the reduced D // and the reduced G that is orthogonal to the reduced D. } case 120: { printState(120); // XXX ++iterc; temp = gredsq * dredsq - dredg * dredg; if (temp <= qred * 1e-4 * qred) { state = 190; break; } temp = Math.sqrt(temp); for (int i = 0; i < n; i++) { if (xbdi.getEntry(i) == ZERO) { s.setEntry(i, (dredg * trialStepPoint.getEntry(i) - dredsq * gnew.getEntry(i)) / temp); } else { s.setEntry(i, ZERO); } } sredg = -temp; // By considering the simple bounds on the variables, calculate an upper // bound on the tangent of half the angle of the alternative iteration, // namely ANGBD, except that, if already a free variable has reached a // bound, there is a branch back to label 100 after fixing that variable. angbd = ONE; iact = -1; for (int i = 0; i < n; i++) { if (xbdi.getEntry(i) == ZERO) { tempa = trustRegionCenterOffset.getEntry(i) + trialStepPoint.getEntry(i) - lowerDifference.getEntry(i); tempb = upperDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i) - trialStepPoint.getEntry(i); if (tempa <= ZERO) { ++nact; xbdi.setEntry(i, MINUS_ONE); state = 100; break; } else if (tempb <= ZERO) { ++nact; xbdi.setEntry(i, ONE); state = 100; break; } // Computing 2nd power double d1 = trialStepPoint.getEntry(i); // Computing 2nd power double d2 = s.getEntry(i); ssq = d1 * d1 + d2 * d2; // Computing 2nd power d1 = trustRegionCenterOffset.getEntry(i) - lowerDifference.getEntry(i); temp = ssq - d1 * d1; if (temp > ZERO) { temp = Math.sqrt(temp) - s.getEntry(i); if (angbd * temp > tempa) { angbd = tempa / temp; iact = i; xsav = MINUS_ONE; } } // Computing 2nd power d1 = upperDifference.getEntry(i) - trustRegionCenterOffset.getEntry(i); temp = ssq - d1 * d1; if (temp > ZERO) { temp = Math.sqrt(temp) + s.getEntry(i); if (angbd * temp > tempb) { angbd = tempb / temp; iact = i; xsav = ONE; } } } } // Calculate HHD and some curvatures for the alternative iteration. state = 210; break; } case 150: { printState(150); // XXX shs = ZERO; dhs = ZERO; dhd = ZERO; for (int i = 0; i < n; i++) { if (xbdi.getEntry(i) == ZERO) { shs += s.getEntry(i) * hs.getEntry(i); dhs += trialStepPoint.getEntry(i) * hs.getEntry(i); dhd += trialStepPoint.getEntry(i) * hred.getEntry(i); } } // Seek the greatest reduction in Q for a range of equally spaced values // of ANGT in [0,ANGBD], where ANGT is the tangent of half the angle of // the alternative iteration. redmax = ZERO; isav = -1; redsav = ZERO; iu = (int) (angbd * 17. + 3.1); for (int i = 0; i < iu; i++) { angt = angbd * i / iu; sth = (angt + angt) / (ONE + angt * angt); temp = shs + angt * (angt * dhd - dhs - dhs); rednew = sth * (angt * dredg - sredg - HALF * sth * temp); if (rednew > redmax) { redmax = rednew; isav = i; rdprev = redsav; } else if (i == isav + 1) { rdnext = rednew; } redsav = rednew; } // Return if the reduction is zero. Otherwise, set the sine and cosine // of the angle of the alternative iteration, and calculate SDEC. if (isav < 0) { state = 190; break; } if (isav < iu) { temp = (rdnext - rdprev) / (redmax + redmax - rdprev - rdnext); angt = angbd * (isav + HALF * temp) / iu; } cth = (ONE - angt * angt) / (ONE + angt * angt); sth = (angt + angt) / (ONE + angt * angt); temp = shs + angt * (angt * dhd - dhs - dhs); sdec = sth * (angt * dredg - sredg - HALF * sth * temp); if (sdec <= ZERO) { state = 190; break; } // Update GNEW, D and HRED. If the angle of the alternative iteration // is restricted by a bound on a free variable, that variable is fixed // at the bound. dredg = ZERO; gredsq = ZERO; for (int i = 0; i < n; i++) { gnew.setEntry(i, gnew.getEntry(i) + (cth - ONE) * hred.getEntry(i) + sth * hs.getEntry(i)); if (xbdi.getEntry(i) == ZERO) { trialStepPoint.setEntry(i, cth * trialStepPoint.getEntry(i) + sth * s.getEntry(i)); dredg += trialStepPoint.getEntry(i) * gnew.getEntry(i); // Computing 2nd power final double d1 = gnew.getEntry(i); gredsq += d1 * d1; } hred.setEntry(i, cth * hred.getEntry(i) + sth * hs.getEntry(i)); } qred += sdec; if (iact >= 0 && isav == iu) { ++nact; xbdi.setEntry(iact, xsav); state = 100; break; } // If SDEC is sufficiently small, then RETURN after setting XNEW to // XOPT+D, giving careful attention to the bounds. if (sdec > qred * .01) { state = 120; break; } } case 190: { printState(190); // XXX dsq = ZERO; for (int i = 0; i < n; i++) { // Computing MAX // Computing MIN final double min = Math.min(trustRegionCenterOffset.getEntry(i) + trialStepPoint.getEntry(i), upperDifference.getEntry(i)); newPoint.setEntry(i, Math.max(min, lowerDifference.getEntry(i))); if (xbdi.getEntry(i) == MINUS_ONE) { newPoint.setEntry(i, lowerDifference.getEntry(i)); } if (xbdi.getEntry(i) == ONE) { newPoint.setEntry(i, upperDifference.getEntry(i)); } trialStepPoint.setEntry(i, newPoint.getEntry(i) - trustRegionCenterOffset.getEntry(i)); // Computing 2nd power final double d1 = trialStepPoint.getEntry(i); dsq += d1 * d1; } return new double[] { dsq, crvmin }; // The following instructions multiply the current S-vector by the second // derivative matrix of the quadratic model, putting the product in HS. // They are reached from three different parts of the software above and // they can be regarded as an external subroutine. } case 210: { printState(210); // XXX int ih = 0; for (int j = 0; j < n; j++) { hs.setEntry(j, ZERO); for (int i = 0; i <= j; i++) { if (i < j) { hs.setEntry(j, hs.getEntry(j) + modelSecondDerivativesValues.getEntry(ih) * s.getEntry(i)); } hs.setEntry(i, hs.getEntry(i) + modelSecondDerivativesValues.getEntry(ih) * s.getEntry(j)); ih++; } } final RealVector tmp = interpolationPoints.operate(s).ebeMultiply(modelSecondDerivativesParameters); for (int k = 0; k < npt; k++) { if (modelSecondDerivativesParameters.getEntry(k) != ZERO) { for (int i = 0; i < n; i++) { hs.setEntry(i, hs.getEntry(i) + tmp.getEntry(k) * interpolationPoints.getEntry(k, i)); } } } if (crvmin != ZERO) { state = 50; break; } if (iterc > itcsav) { state = 150; break; } for (int i = 0; i < n; i++) { hred.setEntry(i, hs.getEntry(i)); } state = 120; break; } default: { throw new MathIllegalStateException(LocalizedFormats.SIMPLE_MESSAGE, "trsbox"); } } } }
From source file:edu.stanford.cfuller.colocalization3d.correction.PositionCorrector.java
/** * Determines the target registration error for a correction by successively leaving out each ImageObject in a set used to make a correction, * calculating a correction from the remaining objects, and assessing the error in correcting the object left out. * //from w w w . jav a 2 s . c o m * @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 The average value of the error over all objects. */ public double determineTRE(java.util.List<ImageObject> imageObjects) { int referenceChannel = this.parameters.getIntValueForKey(REF_CH_PARAM); int channelToCorrect = this.parameters.getIntValueForKey(CORR_CH_PARAM); RealVector treVector = new ArrayRealVector(imageObjects.size(), 0.0); RealVector treXYVector = new ArrayRealVector(imageObjects.size(), 0.0); java.util.Deque<TREThread> startedThreads = new java.util.LinkedList<TREThread>(); int maxThreads = 1; if (this.parameters.hasKey(THREAD_COUNT_PARAM)) { maxThreads = this.parameters.getIntValueForKey(THREAD_COUNT_PARAM); } final int threadWaitTime_ms = 1000; for (int removeIndex = 0; removeIndex < imageObjects.size(); removeIndex++) { if (removeIndex % 10 == 0) { java.util.logging.Logger .getLogger(edu.stanford.cfuller.colocalization3d.Colocalization3DMain.LOGGER_NAME) .finer("calulating TRE: point " + (removeIndex + 1) + " of " + imageObjects.size()); } TREThread nextFit = new TREThread(imageObjects, referenceChannel, channelToCorrect, removeIndex, this); if (startedThreads.size() < maxThreads) { startedThreads.add(nextFit); nextFit.start(); } else { TREThread next = startedThreads.poll(); try { next.join(threadWaitTime_ms); } catch (InterruptedException e) { e.printStackTrace(); } while (next.isAlive()) { startedThreads.add(next); next = startedThreads.poll(); try { next.join(threadWaitTime_ms); } catch (InterruptedException e) { e.printStackTrace(); } } treVector.setEntry(next.getRemoveIndex(), next.getTre()); treXYVector.setEntry(next.getRemoveIndex(), next.getTreXY()); startedThreads.add(nextFit); nextFit.start(); } } java.util.List<Integer> unsuccessful_TRE = new java.util.ArrayList<Integer>(); while (startedThreads.size() > 0) { TREThread next = startedThreads.poll(); try { next.join(); if (next.getSuccess()) { treVector.setEntry(next.getRemoveIndex(), next.getTre()); } else { unsuccessful_TRE.add(next.getRemoveIndex()); } } catch (InterruptedException e) { e.printStackTrace(); } } RealVector treVector_mod = new ArrayRealVector(treVector.getDimension() - unsuccessful_TRE.size()); RealVector treXYVector_mod = new ArrayRealVector(treVector_mod.getDimension()); int c = 0; //unsuccessful TRE calculation results when there is incomplete coverage in the correction dataset for (int i = 0; i < treVector.getDimension(); ++i) { if (!unsuccessful_TRE.contains(i)) { treVector_mod.setEntry(c, treVector.getEntry(i)); treXYVector_mod.setEntry(c, treXYVector.getEntry(i)); ++c; } } treVector = treVector_mod; treXYVector = treXYVector_mod; double tre = treVector.getL1Norm() / treVector.getDimension(); double xy_tre = (treXYVector.getL1Norm() / treXYVector.getDimension()); java.util.logging.Logger.getLogger(edu.stanford.cfuller.colocalization3d.Colocalization3DMain.LOGGER_NAME) .info("TRE: " + tre); java.util.logging.Logger.getLogger(edu.stanford.cfuller.colocalization3d.Colocalization3DMain.LOGGER_NAME) .info("x-y TRE: " + xy_tre); return tre; }
From source file:edu.utexas.cs.tactex.utilityestimation.UtilityEstimatorDefaultForConsumption.java
/** * Core method for estimating utility//from w w w . ja v a2s.c om */ public double estimateUtility(HashMap<TariffSpecification, HashMap<CustomerInfo, Integer>> tariffSubscriptions, HashMap<TariffSpecification, HashMap<CustomerInfo, Double>> predictedCustomerSubscriptions, HashMap<CustomerInfo, HashMap<TariffSpecification, Double>> customer2estimatedTariffCharges, HashMap<CustomerInfo, HashMap<TariffSpecification, ShiftedEnergyData>> customerTariff2ShiftedEnergy, HashMap<CustomerInfo, ArrayRealVector> customer2NonShiftedEnergy, MarketPredictionManager marketPredictionManager, CostCurvesPredictor costCurvesPredictor, int currentTimeslot/*, HashMap<CustomerInfo,ArrayRealVector> customer2estimatedEnergy*/) { // <= for print purposes log.debug( "estimateUtility(): currently assuming competing tariffs are not new and that my subscriptions are not going to change as a result of them"); // accumulate the final results //int expectedRecordLength = 7 * 24; // TODO use Timeservice here correctly double estTariffIncome = 0; double estTariffCosts = 0; double wholesaleCosts = 0; double balancingCosts = 0; double distributionCosts = 0; double withdrawCosts = 0; //double totalConsumption = 0; //double totalProduction = 0; final int predictionRecordLength = BrokerUtils.extractPredictionRecordLength(customerTariff2ShiftedEnergy); RealVector predictedEnergyRecord = new ArrayRealVector(predictionRecordLength); RealVector totalConsumptionEnergyRecord = new ArrayRealVector(predictionRecordLength); RealVector totalProductionEnergyRecord = new ArrayRealVector(predictionRecordLength); // for each customer get his usage prediction for (Entry<TariffSpecification, HashMap<CustomerInfo, Double>> entry : predictedCustomerSubscriptions .entrySet()) { TariffSpecification spec = entry.getKey(); HashMap<CustomerInfo, Double> subscriptions = entry.getValue(); for (Entry<CustomerInfo, Double> ce : subscriptions.entrySet()) { CustomerInfo customerInfo = ce.getKey(); Double subscribedPopulation = ce.getValue(); // Predicted total tariff cash flow. Sign is inverted since // evaluatedTariffs was computed from customers' perspective if (spec.getPowerType().isConsumption()) { estTariffIncome += -customer2estimatedTariffCharges.get(customerInfo).get(spec) * subscribedPopulation; } else if (spec.getPowerType().isProduction()) { estTariffCosts += -customer2estimatedTariffCharges.get(customerInfo).get(spec) * subscribedPopulation; } else { log.warn("Ignoring unknown powertype when computing tariffs income/costs: " + spec.getPowerType()); } // Predicted total energy RealVector energyPrediction = customerTariff2ShiftedEnergy.get(customerInfo).get(spec) .getShiftedEnergy().mapMultiply(subscribedPopulation); predictedEnergyRecord = predictedEnergyRecord.add(energyPrediction); // Predicted balancing cost balancingCosts += 0; log.debug("Ignoring balancing costs - assuming they are 0"); // Predicted withdraw costs (currently assuming everyone will pay). // sign is inverted since withdraw payment is from customer's perspective HashMap<CustomerInfo, Integer> cust2subs = tariffSubscriptions.get(spec); if (cust2subs != null) { Integer currentSubs = cust2subs.get(customerInfo); if (currentSubs != null && subscribedPopulation < currentSubs) { double withdraws = currentSubs - subscribedPopulation; withdrawCosts += -(withdraws * spec.getEarlyWithdrawPayment()); } } // Predicted total consumption and total production if (spec.getPowerType().isConsumption()) { totalConsumptionEnergyRecord = totalConsumptionEnergyRecord.add(energyPrediction); } else if (spec.getPowerType().isProduction()) { totalProductionEnergyRecord = totalProductionEnergyRecord.add(energyPrediction); } else { log.warn("Ignoring unknown powertype when computing distribution costs"); } } } // Predicted distribution costs log.debug("Ignoring balancing orders and curtailment when computing distribution costs"); distributionCosts = 0; double distributionFee = contextManager.getDistributionFee(); for (int i = 0; i < totalConsumptionEnergyRecord.getDimension(); ++i) { double totalTimeslotConsumption = Math.abs(totalConsumptionEnergyRecord.getEntry(i)); double totalTimeslotProduction = Math.abs(totalProductionEnergyRecord.getEntry(i)); distributionCosts += Math.max(totalTimeslotConsumption, totalTimeslotProduction) * distributionFee; } // Predicted wholesale costs (in one of the following two methods:) if (configuratorFactoryService.isUseCostCurves()) { // TODO: might not work for non-fixed rate competitor tariffs - // better to send a mapping of competitor tariffs inside // predictedCustomerSubscriptions, compute shifted predictions for them, // and use these predictions here. // compute energy of customers I don't have RealVector predictedCompetitorsEnergyRecord = new ArrayRealVector(predictionRecordLength); HashMap<CustomerInfo, HashMap<TariffSpecification, Double>> map = BrokerUtils .revertKeyMapping(predictedCustomerSubscriptions); for (CustomerInfo cust : map.keySet()) { double subsToOthers = cust.getPopulation() - BrokerUtils.sumMapValues(map.get(cust)); RealVector customerNonShiftedEnergy = customer2NonShiftedEnergy.get(cust).mapMultiply(subsToOthers); predictedCompetitorsEnergyRecord = predictedCompetitorsEnergyRecord.add(customerNonShiftedEnergy); } for (int i = 0; i < predictedEnergyRecord.getDimension(); ++i) { int futureTimeslot = currentTimeslot + i; double neededKwh = predictedEnergyRecord.getEntry(i); double competitorKwh = predictedCompetitorsEnergyRecord.getEntry(i); double unitCost = costCurvesPredictor.predictUnitCostKwh(currentTimeslot, futureTimeslot, neededKwh, competitorKwh); // NOTE: unitCost is signed (typically negative) wholesaleCosts += (unitCost + costCurvesPredictor.getFudgeFactorKwh(currentTimeslot)) * neededKwh; log.debug("cost-curve prediction: current " + currentTimeslot + " futureTimeslot " + futureTimeslot + " neededKwh " + neededKwh + " neededKwh + competitorKwh " + (neededKwh + competitorKwh) + " unitCost " + unitCost); } } else { // compute wholesale costs using this consumption and market manager marketprices prediction ArrayRealVector estimatedMarketPrices = marketPredictionManager.getPricePerKwhPredictionForAbout7Days(); // sanity check if (predictedEnergyRecord.getDimension() != estimatedMarketPrices.getDimension()) { log.error("Cannot compute utility - prediction periods of market and energy differ);"); return 0; } wholesaleCosts = -predictedEnergyRecord.dotProduct(estimatedMarketPrices); } log.info("estTariffIncome " + estTariffIncome); log.info("estTariffCosts " + estTariffCosts); log.info("wholesaleCosts " + wholesaleCosts); log.info("balancingCosts " + balancingCosts); log.info("distributionCosts " + distributionCosts); //log.info("old distributionCosts " + Math.max(totalProduction, totalConsumption) * contextManager.getDistributionFee()); log.info("withdrawCosts " + withdrawCosts); return estTariffIncome + estTariffCosts + wholesaleCosts + balancingCosts + distributionCosts + withdrawCosts; }