List of usage examples for org.apache.commons.math3.stat.descriptive.moment StandardDeviation StandardDeviation
public StandardDeviation()
From source file:com.dasasian.chok.testutil.loadtest.LoadTestMasterOperation.java
@Override public void nodeOperationsComplete(MasterContext context, List<OperationResult> nodeResults) throws Exception { try {//from ww w .ja va2 s . co m final int queryRate = calculateCurrentQueryRate(); LOG.info("collecting results for iteration " + currentIteration + " and query rate " + queryRate + " after " + (System.currentTimeMillis() - currentIterationStartTime) + " ms ..."); List<LoadTestQueryResult> queryResults = new ArrayList<>(); for (OperationResult operationResult : nodeResults) { if (operationResult == null || operationResult.getUnhandledError() != null) { Exception rootException = null; if (operationResult != null) { //rootException = operationResult.getUnhandledError(); } throw new IllegalStateException( "at least one node operation did not completed properly: " + nodeResults, rootException); } LoadTestNodeOperationResult nodeOperationResult = (LoadTestNodeOperationResult) operationResult; queryResults.addAll(nodeOperationResult.getQueryResults()); } LOG.info("Received " + queryResults.size() + " queries, expected " + queryRate * runTime / 1000); File statisticsFile = new File(resultDir, "load-test-log-" + startTime + ".log"); File resultsFile = new File(resultDir, "load-test-results-" + startTime + ".log"); Writer statisticsWriter = new OutputStreamWriter(new FileOutputStream(statisticsFile, true)); Writer resultWriter = new OutputStreamWriter(new FileOutputStream(resultsFile, true)); if (currentIteration == 0) { // print headers statisticsWriter.append("#queryRate \tnode \tstartTime \tendTime \telapseTime \tquery \n"); resultWriter.append( "#requestedQueryRate \tachievedQueryRate \tfiredQueries \tqueryErrors \tavarageQueryDuration \tstandardDeviation \n"); } try { StorelessUnivariateStatistic timeStandardDeviation = new StandardDeviation(); StorelessUnivariateStatistic timeMean = new Mean(); int errors = 0; for (LoadTestQueryResult result : queryResults) { long elapsedTime = result.getEndTime() > 0 ? result.getEndTime() - result.getStartTime() : -1; statisticsWriter.write(queryRate + "\t" + result.getNodeId() + "\t" + result.getStartTime() + "\t" + result.getEndTime() + "\t" + elapsedTime + "\t" + result.getQuery() + "\n"); if (elapsedTime != -1) { timeStandardDeviation.increment(elapsedTime); timeMean.increment(elapsedTime); } else { ++errors; } } resultWriter.write(queryRate + "\t" + ((double) queryResults.size() / (runTime / 1000)) + "\t" + queryResults.size() + "\t" + errors + "\t" + (int) timeMean.getResult() + "\t" + (int) timeStandardDeviation.getResult() + "\n"); } catch (IOException e) { throw new IllegalStateException("Failed to write statistics data.", e); } try { LOG.info("results written to " + resultsFile.getAbsolutePath()); LOG.info("statistics written to " + statisticsFile.getAbsolutePath()); statisticsWriter.close(); resultWriter.close(); } catch (IOException e) { LOG.warn("Failed to close statistics file."); } if (queryRate + step <= endRate) { currentIteration++; LOG.info("triggering next iteration " + currentIteration); context.getMasterQueue().add(this); } else { LOG.info("finish load test in iteration " + currentIteration + " after " + (System.currentTimeMillis() - startTime) + " ms"); context.getProtocol().removeFlag(getName()); } } catch (Exception e) { context.getProtocol().removeFlag(getName()); } }
From source file:com.musicg.api.DetectionApi.java
protected boolean isPassedStandardDeviation(double[][] spectrogramData) { // normalize the spectrogramData (with all frames in the spectrogram) normalizeSpectrogramData(spectrogramData); // analyst data in this frame // since fftSampleSize==numSamples, there're only one spectrum which is // spectrogramData[last] double[] spectrum = spectrogramData[spectrogramData.length - 1]; // find top most robust frequencies in this frame double[] robustFrequencies = new double[numRobust]; ArrayRankDouble arrayRankDouble = new ArrayRankDouble(); double nthValue = arrayRankDouble.getNthOrderedValue(spectrum, numRobust, false); // end analyst data in this frame int count = 0; for (int i = 0; i < spectrum.length; i++) { if (spectrum[i] >= nthValue) { robustFrequencies[count++] = spectrum[i]; if (count >= numRobust) { break; }//from w w w . ja v a 2s. co m } } // end find top most robust frequencies StandardDeviation standardDeviation = new StandardDeviation(); double sd = standardDeviation.evaluate(robustFrequencies); // range of standard deviation boolean result = (sd >= minStandardDeviation && sd <= maxStandardDeviation); // System.out.println("sd: " + sd + " " + result); return result; }
From source file:gamlss.distributions.ST1.java
/** Calculate and set initial value of sigma. * @param y - vector of values of response variable * @return vector of initial values of sigma *//* w ww . jav a 2 s . c om*/ private ArrayRealVector setSigmaInitial(final ArrayRealVector y) { //sigma.initial = expression(sigma <- rep(sd(y)/4, length(y))), tempV = new ArrayRealVector(y.getDimension()); final double out = new StandardDeviation().evaluate(y.getDataRef()); tempV.set(out / 4); return tempV; }
From source file:gedi.util.math.stat.distributions.NormalMixtureDistribution.java
public static NormalMixtureDistribution init(final double[] data, final int numComponents) throws NotStrictlyPositiveException, DimensionMismatchException { if (numComponents == 1) return new NormalMixtureDistribution(new NormalDistribution[] { new NormalDistribution(new Mean().evaluate(data), new StandardDeviation().evaluate(data)) }, new double[] { 1 }); if (data.length < 2) { throw new NotStrictlyPositiveException(data.length); }//from w w w . j a v a2 s .c o m if (numComponents < 1) { throw new NumberIsTooSmallException(numComponents, 2, true); } if (numComponents > data.length) { throw new NumberIsTooLargeException(numComponents, data.length, true); } final int numRows = data.length; double[] sortedData = data.clone(); Arrays.sort(sortedData); // components of mixture model to be created double[] mixing = new double[numComponents]; NormalDistribution[] comp = new NormalDistribution[numComponents]; // create a component based on data in each bin for (int k = 0; k < numComponents; k++) { // minimum index (inclusive) from sorted data for this bin final int minIndex = (k * numRows) / numComponents; // maximum index (exclusive) from sorted data for this bin final int maxIndex = Math.min(numRows, ((k + 1) * numRows) / numComponents); double m = new Mean().evaluate(sortedData, minIndex, maxIndex - minIndex); double sd = new StandardDeviation().evaluate(sortedData, minIndex, maxIndex - minIndex); mixing[k] = 1d / numComponents; comp[k] = new NormalDistribution(m, sd); } return new NormalMixtureDistribution(comp, mixing); }
From source file:com.itemanalysis.jmetrik.graph.nicc.NonparametricCurveAnalysis.java
private void initializeGridPoints() throws SQLException { Statement stmt = null;/*from www .j a va2s .c o m*/ ResultSet rs = null; //connect to db try { Table sqlTable = new Table(tableName.getNameForDatabase()); SelectQuery select = new SelectQuery(); select.addColumn(sqlTable, regressorVariable.getName().nameForDatabase()); stmt = conn.createStatement(ResultSet.TYPE_SCROLL_INSENSITIVE, ResultSet.CONCUR_READ_ONLY); rs = stmt.executeQuery(select.toString()); Min min = new Min(); Max max = new Max(); Mean mean = new Mean(); StandardDeviation sd = new StandardDeviation(); double value = 0.0; while (rs.next()) { value = rs.getDouble(regressorVariable.getName().nameForDatabase()); if (!rs.wasNull()) { min.increment(value); max.increment(value); mean.increment(value); sd.increment(value); } updateProgress(); } rs.close(); stmt.close(); //evaluation points double sdv = sd.getResult(); double mn = mean.getResult(); double lower = mn - 2.5 * sdv; double upper = mn + 2.5 * sdv; bwAdjustment *= sdv; bandwidth = new NonparametricIccBandwidth(sampleSize, bwAdjustment); gridPoints = command.getFreeOption("gridpoints").getInteger(); // uniformDistributionApproximation = new UniformDistributionApproximation( // min.getResult(), max.getResult(), gridPoints); uniformDistributionApproximation = new UniformDistributionApproximation(lower, upper, gridPoints); } catch (SQLException ex) { throw ex; } finally { if (rs != null) rs.close(); if (stmt != null) stmt.close(); } }
From source file:de.biomedical_imaging.ij.nanotrackj.Track.java
/** * //w w w .jav a2 s. c o m * @param driftx Drift in x direction (in pixels) * @param drifty Drift in y direction (in pixels) * @param tau Timelag * @return [0] = The mean squared displacement, [1] = Number of data points */ public double[] getMeanSquareDisplacementSD(double driftx, double drifty, int tau) { double msd = 0; StandardDeviation sd = new StandardDeviation(); int N = 0; for (int i = tau; i < this.size(); ++i) { msd = Math.pow(this.get(i - tau).getX() - this.get(i).getX() - tau * driftx, 2) + Math.pow(this.get(i - tau).getY() - this.get(i).getY() - tau * drifty, 2); sd.increment(msd); ++N; } double[] result = new double[2]; result[0] = sd.getResult(); result[1] = N; return result; }
From source file:net.sf.javaml.clustering.AQBC.java
/** * Normalizes the data to mean 0 and standard deviation 1. This method * discards all instances that cannot be normalized, i.e. they have the same * value for all attributes.// www .j ava 2 s.co m * * @param data * @return */ private Vector<TaggedInstance> normalize(Dataset data) { Vector<TaggedInstance> out = new Vector<TaggedInstance>(); for (int i = 0; i < data.size(); i++) { Double[] old = data.instance(i).values().toArray(new Double[0]); double[] conv = new double[old.length]; for (int j = 0; j < old.length; j++) { conv[j] = old[j]; } Mean m = new Mean(); double MU = m.evaluate(conv); // System.out.println("MU = "+MU); StandardDeviation std = new StandardDeviation(); double SIGM = std.evaluate(conv, MU); // System.out.println("SIGM = "+SIGM); if (!MathUtils.eq(SIGM, 0)) { double[] val = new double[old.length]; for (int j = 0; j < old.length; j++) { val[j] = (float) ((old[j] - MU) / SIGM); } // System.out.println("VAL "+i+" = "+Arrays.toString(val)); out.add(new TaggedInstance(new DenseInstance(val, data.instance(i).classValue()), i)); } } // System.out.println("FIRST = "+out.get(0)); return out; }
From source file:de.biomedical_imaging.traj.math.TrajectorySplineFitLegacy.java
/** * Calculates a spline to a trajectory. Attention: The spline is fitted through a rotated version of the trajectory. * To fit the spline the trajectory is rotated into its main direction. You can access this rotated trajectory by * {@link #getRotatedTrajectory() getRotatedTrajectory}. * @return The fitted spline//from w ww.ja v a 2 s . c o m */ public PolynomialSplineFunction calculateSpline() { /* * 1.Calculate the minimum bounding rectangle */ ArrayList<Point2D.Double> points = new ArrayList<Point2D.Double>(); for (int i = 0; i < t.size(); i++) { Point2D.Double p = new Point2D.Double(); p.setLocation(t.get(i).x, t.get(i).y); points.add(p); } Point2D.Double[] rect = null; try { rect = RotatingCalipers.getMinimumBoundingRectangle(points); } catch (IllegalArgumentException e) { } catch (EmptyStackException e) { } /* * 1.1 Rotate that the major axis is parallel with the xaxis */ Point2D.Double majorDirection = null; Point2D.Double p1 = rect[2]; //top left Point2D.Double p2 = p1.distance(rect[1]) > p1.distance(rect[3]) ? rect[1] : rect[3]; //Point to long side Point2D.Double p3 = p1.distance(rect[1]) > p1.distance(rect[3]) ? rect[3] : rect[1]; //Point to short side majorDirection = new Point2D.Double(p2.x - p1.x, p2.y - p1.y); double width = p1.distance(p2); double inRad = -1 * Math.atan2(majorDirection.y, majorDirection.x); boolean doTransform = (Math.abs(Math.abs(inRad) - Math.PI) > 0.001); if (doTransform) { angleRotated = inRad; for (int i = 0; i < t.size(); i++) { double x = t.get(i).x; double y = t.get(i).y; double newX = x * Math.cos(inRad) - y * Math.sin(inRad); double newY = x * Math.sin(inRad) + y * Math.cos(inRad); rotatedTrajectory.add(newX, newY, 0); points.get(i).setLocation(newX, newY); } for (int i = 0; i < rect.length; i++) { rect[i].setLocation(rect[i].x * Math.cos(inRad) - rect[i].y * Math.sin(inRad), rect[i].x * Math.sin(inRad) + rect[i].y * Math.cos(inRad)); } p1 = rect[2]; //top left p2 = p1.distance(rect[1]) > p1.distance(rect[3]) ? rect[1] : rect[3]; //Point to long side p3 = p1.distance(rect[1]) > p1.distance(rect[3]) ? rect[3] : rect[1]; //Point to short side } else { angleRotated = 0; rotatedTrajectory = t; } /* * 2. Divide the rectangle in n equal segments * 2.1 Calculate line in main direction * 2.2 Project the points in onto this line * 2.3 Calculate the distance between the start of the line and the projected point * 2.4 Assign point to segment according to distance of (2.3) */ List<List<Point2D.Double>> pointsInSegments = null; boolean allSegmentsContainingAtLeastTwoPoints = true; do { allSegmentsContainingAtLeastTwoPoints = true; double segmentWidth = p1.distance(p2) / nSegments; pointsInSegments = new ArrayList<List<Point2D.Double>>(nSegments); for (int i = 0; i < nSegments; i++) { pointsInSegments.add(new ArrayList<Point2D.Double>()); } for (int i = 0; i < points.size(); i++) { Point2D.Double projPoint = projectPointToLine(p1, p2, points.get(i)); int index = (int) (p1.distance(projPoint) / segmentWidth); if (index > (nSegments - 1)) { index = (nSegments - 1); } pointsInSegments.get(index).add(points.get(i)); } for (int i = 0; i < pointsInSegments.size(); i++) { if (pointsInSegments.get(i).size() < 2) { if (nSegments > 2) { nSegments--; i = pointsInSegments.size(); allSegmentsContainingAtLeastTwoPoints = false; } } } } while (allSegmentsContainingAtLeastTwoPoints == false); /* * 3. Calculate the mean standard deviation over each segment: <s> */ Point2D.Double eMajorP1 = new Point2D.Double(p1.x - (p3.x - p1.x) / 2.0, p1.y - (p3.y - p1.y) / 2.0); Point2D.Double eMajorP2 = new Point2D.Double(p2.x - (p3.x - p1.x) / 2.0, p2.y - (p3.y - p1.y) / 2.0); double sumMean = 0; int Nsum = 0; for (int i = 0; i < nSegments; i++) { StandardDeviation sd = new StandardDeviation(); double[] distances = new double[pointsInSegments.get(i).size()]; for (int j = 0; j < pointsInSegments.get(i).size(); j++) { int factor = 1; if (isLeft(eMajorP1, eMajorP2, pointsInSegments.get(i).get(j))) { factor = -1; } distances[j] = factor * distancePointLine(eMajorP1, eMajorP2, pointsInSegments.get(i).get(j)); } if (distances.length > 0) { sd.setData(distances); sumMean += sd.evaluate(); Nsum++; } } double s = sumMean / Nsum; if (s < 0.000000000001) { s = width / nSegments; } /* * 4. Build a kd-tree */ KDTree<Point2D.Double> kdtree = new KDTree<Point2D.Double>(2); for (int i = 0; i < points.size(); i++) { try { //To ensure that all points have a different key, add small random number kdtree.insert(new double[] { points.get(i).x, points.get(i).y }, points.get(i)); } catch (KeySizeException e) { e.printStackTrace(); } catch (KeyDuplicateException e) { //Do nothing! It is not important } } /* * 5. Using the first point f in trajectory and calculate the center of mass * of all points around f (radius: 3*<s>)) */ List<Point2D.Double> near = null; Point2D.Double first = minDistancePointToLine(p1, p3, points); double r1 = 3 * s; try { near = kdtree.nearestEuclidean(new double[] { first.x, first.y }, r1); } catch (KeySizeException e) { e.printStackTrace(); } double cx = 0; double cy = 0; for (int i = 0; i < near.size(); i++) { cx += near.get(i).x; cy += near.get(i).y; } cx /= near.size(); cy /= near.size(); splineSupportPoints = new ArrayList<Point2D.Double>(); splineSupportPoints.add(new Point2D.Double(cx, cy)); /* * 6. The second point is determined by finding the center-of-mass of particles in the p/2 radian * section of an annulus, r1 < r < 2r1, that is directed toward the angle with the highest number * of particles within p/2 radians. * 7. This second point is then used as the center of the annulus for choosing the third point, and the process is repeated (6. & 7.). */ /* * 6.1 Find all points in the annolous */ /* * 6.2 Write each point in a coordinate system centered at the center of the sphere, calculate direction and * check if it in the allowed bounds */ int nCircleSegments = 100; double deltaRad = 2 * Math.PI / nCircleSegments; boolean stop = false; int minN = 7; double tempr1 = r1; double allowedDeltaDirection = 0.5 * Math.PI; while (stop == false) { List<Point2D.Double> nearestr1 = null; List<Point2D.Double> nearest2xr1 = null; try { nearestr1 = kdtree .nearestEuclidean(new double[] { splineSupportPoints.get(splineSupportPoints.size() - 1).x, splineSupportPoints.get(splineSupportPoints.size() - 1).y }, tempr1); nearest2xr1 = kdtree .nearestEuclidean(new double[] { splineSupportPoints.get(splineSupportPoints.size() - 1).x, splineSupportPoints.get(splineSupportPoints.size() - 1).y }, 2 * tempr1); } catch (KeySizeException e) { // TODO Auto-generated catch block e.printStackTrace(); } nearest2xr1.removeAll(nearestr1); double lThreshRad = 0; double hThreshRad = Math.PI / 2; double stopThresh = 2 * Math.PI; if (splineSupportPoints.size() > 1) { double directionInRad = Math.atan2( splineSupportPoints.get(splineSupportPoints.size() - 1).y - splineSupportPoints.get(splineSupportPoints.size() - 2).y, splineSupportPoints.get(splineSupportPoints.size() - 1).x - splineSupportPoints.get(splineSupportPoints.size() - 2).x) + Math.PI; lThreshRad = directionInRad - allowedDeltaDirection / 2 - Math.PI / 4; if (lThreshRad < 0) { lThreshRad = 2 * Math.PI + lThreshRad; } if (lThreshRad > 2 * Math.PI) { lThreshRad = lThreshRad - 2 * Math.PI; } hThreshRad = directionInRad + allowedDeltaDirection / 2 + Math.PI / 4; if (hThreshRad < 0) { hThreshRad = 2 * Math.PI + hThreshRad; } if (hThreshRad > 2 * Math.PI) { hThreshRad = hThreshRad - 2 * Math.PI; } stopThresh = directionInRad + allowedDeltaDirection / 2 - Math.PI / 4; if (stopThresh > 2 * Math.PI) { stopThresh = stopThresh - 2 * Math.PI; } } double newCx = 0; double newCy = 0; int newCN = 0; int candN = 0; //Find center with highest density of points double lastDist = 0; double newDist = 0; do { lastDist = Math.min(Math.abs(lThreshRad - stopThresh), 2 * Math.PI - Math.abs(lThreshRad - stopThresh)); candN = 0; double candCx = 0; double candCy = 0; for (int i = 0; i < nearest2xr1.size(); i++) { Point2D.Double centerOfCircle = splineSupportPoints.get(splineSupportPoints.size() - 1); Vector2d relativeToCircle = new Vector2d(nearest2xr1.get(i).x - centerOfCircle.x, nearest2xr1.get(i).y - centerOfCircle.y); relativeToCircle.normalize(); double angleInRadians = Math.atan2(relativeToCircle.y, relativeToCircle.x) + Math.PI; if (lThreshRad < hThreshRad) { if (angleInRadians > lThreshRad && angleInRadians < hThreshRad) { candCx += nearest2xr1.get(i).x; candCy += nearest2xr1.get(i).y; candN++; } } else { if (angleInRadians > lThreshRad || angleInRadians < hThreshRad) { candCx += nearest2xr1.get(i).x; candCy += nearest2xr1.get(i).y; candN++; } } } if (candN > 0 && candN > newCN) { candCx /= candN; candCy /= candN; newCx = candCx; newCy = candCy; newCN = candN; } lThreshRad += deltaRad; hThreshRad += deltaRad; if (lThreshRad > 2 * Math.PI) { lThreshRad = lThreshRad - 2 * Math.PI; } if (hThreshRad > 2 * Math.PI) { hThreshRad = hThreshRad - 2 * Math.PI; } newDist = Math.min(Math.abs(lThreshRad - stopThresh), 2 * Math.PI - Math.abs(lThreshRad - stopThresh)); } while ((newDist - lastDist) > 0); //Check if the new center is valid if (splineSupportPoints.size() > 1) { double currentDirectionInRad = Math.atan2( splineSupportPoints.get(splineSupportPoints.size() - 1).y - splineSupportPoints.get(splineSupportPoints.size() - 2).y, splineSupportPoints.get(splineSupportPoints.size() - 1).x - splineSupportPoints.get(splineSupportPoints.size() - 2).x) + Math.PI; double candDirectionInRad = Math.atan2( newCy - splineSupportPoints.get(splineSupportPoints.size() - 1).y, newCx - splineSupportPoints.get(splineSupportPoints.size() - 1).x) + Math.PI; double dDir = Math.max(currentDirectionInRad, candDirectionInRad) - Math.min(currentDirectionInRad, candDirectionInRad); if (dDir > 2 * Math.PI) { dDir = 2 * Math.PI - dDir; } if (dDir > allowedDeltaDirection) { stop = true; } } boolean enoughPoints = (newCN < minN); boolean isNormalRadius = Math.abs(tempr1 - r1) < Math.pow(10, -18); boolean isExtendedRadius = Math.abs(tempr1 - 3 * r1) < Math.pow(10, -18); if (enoughPoints && isNormalRadius) { //Not enough points, extend search radius tempr1 = 3 * r1; } else if (enoughPoints && isExtendedRadius) { //Despite radius extension: Not enough points! stop = true; } else if (stop == false) { splineSupportPoints.add(new Point2D.Double(newCx, newCy)); tempr1 = r1; } } //Sort Collections.sort(splineSupportPoints, new Comparator<Point2D.Double>() { public int compare(Point2D.Double o1, Point2D.Double o2) { if (o1.x < o2.x) { return -1; } if (o1.x > o2.x) { return 1; } return 0; } }); //Add endpoints if (splineSupportPoints.size() > 1) { Vector2d start = new Vector2d(splineSupportPoints.get(0).x - splineSupportPoints.get(1).x, splineSupportPoints.get(0).y - splineSupportPoints.get(1).y); start.normalize(); start.scale(r1 * 8); splineSupportPoints.add(0, new Point2D.Double(splineSupportPoints.get(0).x + start.x, splineSupportPoints.get(0).y + start.y)); Vector2d end = new Vector2d( splineSupportPoints.get(splineSupportPoints.size() - 1).x - splineSupportPoints.get(splineSupportPoints.size() - 2).x, splineSupportPoints.get(splineSupportPoints.size() - 1).y - splineSupportPoints.get(splineSupportPoints.size() - 2).y); end.normalize(); end.scale(r1 * 6); splineSupportPoints .add(new Point2D.Double(splineSupportPoints.get(splineSupportPoints.size() - 1).x + end.x, splineSupportPoints.get(splineSupportPoints.size() - 1).y + end.y)); } else { Vector2d majordir = new Vector2d(-1, 0); majordir.normalize(); majordir.scale(r1 * 8); splineSupportPoints.add(0, new Point2D.Double(splineSupportPoints.get(0).x + majordir.x, splineSupportPoints.get(0).y + majordir.y)); majordir.scale(-1); splineSupportPoints .add(new Point2D.Double(splineSupportPoints.get(splineSupportPoints.size() - 1).x + majordir.x, splineSupportPoints.get(splineSupportPoints.size() - 1).y + majordir.y)); } //Interpolate spline double[] supX = new double[splineSupportPoints.size()]; double[] supY = new double[splineSupportPoints.size()]; for (int i = 0; i < splineSupportPoints.size(); i++) { supX[i] = splineSupportPoints.get(i).x; supY[i] = splineSupportPoints.get(i).y; } SplineInterpolator sIinter = new SplineInterpolator(); spline = sIinter.interpolate(supX, supY); return spline; }
From source file:de.biomedical_imaging.traj.math.TrajectorySplineFit.java
/** * Calculates a spline to a trajectory. Attention: The spline is fitted through a rotated version of the trajectory. * To fit the spline the trajectory is rotated into its main direction. You can access this rotated trajectory by * {@link #getRotatedTrajectory() getRotatedTrajectory}. * @return The fitted spline//from www.j ava 2 s . c om */ public PolynomialSplineFunction calculateSpline() { successfull = false; /* * 1.Calculate the minimum bounding rectangle */ ArrayList<Point2D.Double> points = new ArrayList<Point2D.Double>(); for (int i = 0; i < t.size(); i++) { Point2D.Double p = new Point2D.Double(); p.setLocation(t.get(i).x, t.get(i).y); points.add(p); } /* * 1.1 Rotate that the major axis is parallel with the xaxis */ Array2DRowRealMatrix gyr = RadiusGyrationTensor2D.getRadiusOfGyrationTensor(t); EigenDecomposition eigdec = new EigenDecomposition(gyr); double inRad = -1 * Math.atan2(eigdec.getEigenvector(0).getEntry(1), eigdec.getEigenvector(0).getEntry(0)); boolean doTransform = (Math.abs(Math.abs(inRad) - Math.PI) > 0.001); if (doTransform) { angleRotated = inRad; for (int i = 0; i < t.size(); i++) { double x = t.get(i).x; double y = t.get(i).y; double newX = x * Math.cos(inRad) - y * Math.sin(inRad); double newY = x * Math.sin(inRad) + y * Math.cos(inRad); rotatedTrajectory.add(newX, newY, 0); points.get(i).setLocation(newX, newY); } //for(int i = 0; i < rect.length; i++){ // rect[i].setLocation(rect[i].x*Math.cos(inRad)-rect[i].y*Math.sin(inRad), rect[i].x*Math.sin(inRad)+rect[i].y*Math.cos(inRad)); //} } else { angleRotated = 0; rotatedTrajectory = t; } /* * 2. Divide the rectangle in n equal segments * 2.1 Calculate line in main direction * 2.2 Project the points in onto this line * 2.3 Calculate the distance between the start of the line and the projected point * 2.4 Assign point to segment according to distance of (2.3) */ List<List<Point2D.Double>> pointsInSegments = null; boolean allSegmentsContainingAtLeastTwoPoints = true; int indexSmallestX = 0; double segmentWidth = 0; do { double smallestX = Double.MAX_VALUE; double largestX = Double.MIN_VALUE; for (int i = 0; i < points.size(); i++) { if (points.get(i).x < smallestX) { smallestX = points.get(i).x; indexSmallestX = i; } if (points.get(i).x > largestX) { largestX = points.get(i).x; } } allSegmentsContainingAtLeastTwoPoints = true; segmentWidth = (largestX - smallestX) / nSegments; pointsInSegments = new ArrayList<List<Point2D.Double>>(nSegments); for (int i = 0; i < nSegments; i++) { pointsInSegments.add(new ArrayList<Point2D.Double>()); } for (int i = 0; i < points.size(); i++) { int index = (int) Math.abs((points.get(i).x / segmentWidth)); if (index > (nSegments - 1)) { index = (nSegments - 1); } pointsInSegments.get(index).add(points.get(i)); } for (int i = 0; i < pointsInSegments.size(); i++) { if (pointsInSegments.get(i).size() < 2) { if (nSegments > 2) { nSegments--; i = pointsInSegments.size(); allSegmentsContainingAtLeastTwoPoints = false; } } } } while (allSegmentsContainingAtLeastTwoPoints == false); /* * 3. Calculate the mean standard deviation over each segment: <s> */ //Point2D.Double eMajorP1 = new Point2D.Double(p1.x - (p3.x-p1.x)/2.0,p1.y - (p3.y-p1.y)/2.0); // Point2D.Double eMajorP2 = new Point2D.Double(p2.x - (p3.x-p1.x)/2.0,p2.y - (p3.y-p1.y)/2.0); double sumSDOrthogonal = 0; int Nsum = 0; CenterOfGravityFeature cogf = new CenterOfGravityFeature(rotatedTrajectory); Point2D.Double cog = new Point2D.Double(cogf.evaluate()[0], cogf.evaluate()[1]); Point2D.Double eMajorP1 = cog; Point2D.Double eMajorP2 = new Point2D.Double(cog.getX() + 1, cog.getY()); for (int i = 0; i < nSegments; i++) { StandardDeviation sd = new StandardDeviation(); double[] distances = new double[pointsInSegments.get(i).size()]; for (int j = 0; j < pointsInSegments.get(i).size(); j++) { int factor = 1; if (isLeft(eMajorP1, eMajorP2, pointsInSegments.get(i).get(j))) { factor = -1; } distances[j] = factor * distancePointLine(eMajorP1, eMajorP2, pointsInSegments.get(i).get(j)); } if (distances.length > 0) { sd.setData(distances); sumSDOrthogonal += sd.evaluate(); Nsum++; } } double s = sumSDOrthogonal / Nsum; if (segmentWidth < Math.pow(10, -15)) { return null; } if (s < Math.pow(10, -15)) { //If standard deviation is zero, replace it with the half of the segment width s = segmentWidth / 2; } //rotatedTrajectory.showTrajectory("rot"); /* * 4. Build a kd-tree */ KDTree<Point2D.Double> kdtree = new KDTree<Point2D.Double>(2); for (int i = 0; i < points.size(); i++) { try { //To ensure that all points have a different key, add small random number kdtree.insert(new double[] { points.get(i).x, points.get(i).y }, points.get(i)); } catch (KeySizeException e) { e.printStackTrace(); } catch (KeyDuplicateException e) { //Do nothing! It is not important } } /* * 5. Using the first point f in trajectory and calculate the center of mass * of all points around f (radius: 3*<s>)) */ List<Point2D.Double> near = null; Point2D.Double first = points.get(indexSmallestX);//minDistancePointToLine(p1, p3, points); double r1 = 3 * s; try { near = kdtree.nearestEuclidean(new double[] { first.x, first.y }, r1); } catch (KeySizeException e) { e.printStackTrace(); } double cx = 0; double cy = 0; for (int i = 0; i < near.size(); i++) { cx += near.get(i).x; cy += near.get(i).y; } cx /= near.size(); cy /= near.size(); splineSupportPoints = new ArrayList<Point2D.Double>(); splineSupportPoints.add(new Point2D.Double(cx, cy)); /* * 6. The second point is determined by finding the center-of-mass of particles in the p/2 radian * section of an annulus, r1 < r < 2r1, that is directed toward the angle with the highest number * of particles within p/2 radians. * 7. This second point is then used as the center of the annulus for choosing the third point, and the process is repeated (6. & 7.). */ /* * 6.1 Find all points in the annolous */ /* * 6.2 Write each point in a coordinate system centered at the center of the sphere, calculate direction and * check if it in the allowed bounds */ int nCircleSegments = 100; double deltaRad = 2 * Math.PI / nCircleSegments; boolean stop = false; int minN = 7; double tempr1 = r1; double allowedDeltaDirection = 0.5 * Math.PI; while (stop == false) { List<Point2D.Double> nearestr1 = null; List<Point2D.Double> nearest2xr1 = null; try { nearestr1 = kdtree .nearestEuclidean(new double[] { splineSupportPoints.get(splineSupportPoints.size() - 1).x, splineSupportPoints.get(splineSupportPoints.size() - 1).y }, tempr1); nearest2xr1 = kdtree .nearestEuclidean(new double[] { splineSupportPoints.get(splineSupportPoints.size() - 1).x, splineSupportPoints.get(splineSupportPoints.size() - 1).y }, 2 * tempr1); } catch (KeySizeException e) { // TODO Auto-generated catch block e.printStackTrace(); } nearest2xr1.removeAll(nearestr1); double lThreshRad = 0; double hThreshRad = Math.PI / 2; double stopThresh = 2 * Math.PI; if (splineSupportPoints.size() > 1) { double directionInRad = Math.atan2( splineSupportPoints.get(splineSupportPoints.size() - 1).y - splineSupportPoints.get(splineSupportPoints.size() - 2).y, splineSupportPoints.get(splineSupportPoints.size() - 1).x - splineSupportPoints.get(splineSupportPoints.size() - 2).x) + Math.PI; lThreshRad = directionInRad - allowedDeltaDirection / 2 - Math.PI / 4; if (lThreshRad < 0) { lThreshRad = 2 * Math.PI + lThreshRad; } if (lThreshRad > 2 * Math.PI) { lThreshRad = lThreshRad - 2 * Math.PI; } hThreshRad = directionInRad + allowedDeltaDirection / 2 + Math.PI / 4; if (hThreshRad < 0) { hThreshRad = 2 * Math.PI + hThreshRad; } if (hThreshRad > 2 * Math.PI) { hThreshRad = hThreshRad - 2 * Math.PI; } stopThresh = directionInRad + allowedDeltaDirection / 2 - Math.PI / 4; if (stopThresh > 2 * Math.PI) { stopThresh = stopThresh - 2 * Math.PI; } } double newCx = 0; double newCy = 0; int newCN = 0; int candN = 0; //Find center with highest density of points double lastDist = 0; double newDist = 0; do { lastDist = Math.min(Math.abs(lThreshRad - stopThresh), 2 * Math.PI - Math.abs(lThreshRad - stopThresh)); candN = 0; double candCx = 0; double candCy = 0; for (int i = 0; i < nearest2xr1.size(); i++) { Point2D.Double centerOfCircle = splineSupportPoints.get(splineSupportPoints.size() - 1); Vector2d relativeToCircle = new Vector2d(nearest2xr1.get(i).x - centerOfCircle.x, nearest2xr1.get(i).y - centerOfCircle.y); relativeToCircle.normalize(); double angleInRadians = Math.atan2(relativeToCircle.y, relativeToCircle.x) + Math.PI; if (lThreshRad < hThreshRad) { if (angleInRadians > lThreshRad && angleInRadians < hThreshRad) { candCx += nearest2xr1.get(i).x; candCy += nearest2xr1.get(i).y; candN++; } } else { if (angleInRadians > lThreshRad || angleInRadians < hThreshRad) { candCx += nearest2xr1.get(i).x; candCy += nearest2xr1.get(i).y; candN++; } } } if (candN > 0 && candN > newCN) { candCx /= candN; candCy /= candN; newCx = candCx; newCy = candCy; newCN = candN; } lThreshRad += deltaRad; hThreshRad += deltaRad; if (lThreshRad > 2 * Math.PI) { lThreshRad = lThreshRad - 2 * Math.PI; } if (hThreshRad > 2 * Math.PI) { hThreshRad = hThreshRad - 2 * Math.PI; } newDist = Math.min(Math.abs(lThreshRad - stopThresh), 2 * Math.PI - Math.abs(lThreshRad - stopThresh)); } while ((newDist - lastDist) > 0); //Check if the new center is valid if (splineSupportPoints.size() > 1) { double currentDirectionInRad = Math.atan2( splineSupportPoints.get(splineSupportPoints.size() - 1).y - splineSupportPoints.get(splineSupportPoints.size() - 2).y, splineSupportPoints.get(splineSupportPoints.size() - 1).x - splineSupportPoints.get(splineSupportPoints.size() - 2).x) + Math.PI; double candDirectionInRad = Math.atan2( newCy - splineSupportPoints.get(splineSupportPoints.size() - 1).y, newCx - splineSupportPoints.get(splineSupportPoints.size() - 1).x) + Math.PI; double dDir = Math.max(currentDirectionInRad, candDirectionInRad) - Math.min(currentDirectionInRad, candDirectionInRad); if (dDir > 2 * Math.PI) { dDir = 2 * Math.PI - dDir; } if (dDir > allowedDeltaDirection) { stop = true; } } boolean enoughPoints = (newCN < minN); boolean isNormalRadius = Math.abs(tempr1 - r1) < Math.pow(10, -18); boolean isExtendedRadius = Math.abs(tempr1 - 3 * r1) < Math.pow(10, -18); if (enoughPoints && isNormalRadius) { //Not enough points, extend search radius tempr1 = 3 * r1; } else if (enoughPoints && isExtendedRadius) { //Despite radius extension: Not enough points! stop = true; } else if (stop == false) { splineSupportPoints.add(new Point2D.Double(newCx, newCy)); tempr1 = r1; } } //Sort Collections.sort(splineSupportPoints, new Comparator<Point2D.Double>() { public int compare(Point2D.Double o1, Point2D.Double o2) { if (o1.x < o2.x) { return -1; } if (o1.x > o2.x) { return 1; } return 0; } }); //Add endpoints if (splineSupportPoints.size() > 1) { Vector2d start = new Vector2d(splineSupportPoints.get(0).x - splineSupportPoints.get(1).x, splineSupportPoints.get(0).y - splineSupportPoints.get(1).y); start.normalize(); start.scale(r1 * 8); splineSupportPoints.add(0, new Point2D.Double(splineSupportPoints.get(0).x + start.x, splineSupportPoints.get(0).y + start.y)); Vector2d end = new Vector2d( splineSupportPoints.get(splineSupportPoints.size() - 1).x - splineSupportPoints.get(splineSupportPoints.size() - 2).x, splineSupportPoints.get(splineSupportPoints.size() - 1).y - splineSupportPoints.get(splineSupportPoints.size() - 2).y); end.normalize(); end.scale(r1 * 6); splineSupportPoints .add(new Point2D.Double(splineSupportPoints.get(splineSupportPoints.size() - 1).x + end.x, splineSupportPoints.get(splineSupportPoints.size() - 1).y + end.y)); } else { Vector2d majordir = new Vector2d(-1, 0); majordir.normalize(); majordir.scale(r1 * 8); splineSupportPoints.add(0, new Point2D.Double(splineSupportPoints.get(0).x + majordir.x, splineSupportPoints.get(0).y + majordir.y)); majordir.scale(-1); splineSupportPoints .add(new Point2D.Double(splineSupportPoints.get(splineSupportPoints.size() - 1).x + majordir.x, splineSupportPoints.get(splineSupportPoints.size() - 1).y + majordir.y)); } //Interpolate spline double[] supX = new double[splineSupportPoints.size()]; double[] supY = new double[splineSupportPoints.size()]; for (int i = 0; i < splineSupportPoints.size(); i++) { supX[i] = splineSupportPoints.get(i).x; supY[i] = splineSupportPoints.get(i).y; } SplineInterpolator sIinter = new SplineInterpolator(); spline = sIinter.interpolate(supX, supY); successfull = true; return spline; }
From source file:com.github.rinde.rinsim.scenario.measure.MetricsTest.java
@Test public void dynamismTest2() { final double[] ordersPerHour = { 15d };// , 20d, 50d, 100d, 1000d }; final StandardDeviation sd = new StandardDeviation(); final RandomGenerator rng = new MersenneTwister(123L); final List<Times> times = newArrayList(); // for (int i = 0; i < 10; i++) { // times.add(generateTimes(rng)); // }//from w ww . java 2 s.co m times.add(asTimes(1000, 250L, 500L, 750L)); times.add(asTimes(1000, 100L, 500L, 750L)); times.add(asTimes(1000, 100L, 200L, 300L, 400L, 500L, 600L, 700L, 800L, 900L)); times.add(asTimes(1000, 100L, 200L, 300L, 399L, 500L, 600L, 700L, 800L, 900L)); times.add(asTimes(1000, 10L, 150L, 250L, 350L, 450L, 550L, 650L, 750L, 850L, 950L)); times.add(asTimes(1000, 50L, 150L, 250L, 350L, 450L, 551L, 650L, 750L, 850L, 950L)); times.add(asTimes(1000, 250L, 500L, 750L)); times.add(asTimes(1000, 0L, 50L, 55L, 57L, 59L, 60L, 100L, 150L, 750L)); times.add(asTimes(1000, 5L, 5L, 5L, 5L)); times.add(asTimes(1000, 4L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L)); times.add(asTimes(1000, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L)); times.add(asTimes(1000, 0L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 999L)); times.add(asTimes(1000, 500L, 500L, 500L, 500L)); times.add(asTimes(1000, 5L, 5L, 5L, 5L, 400L, 410L, 430L, 440L, 800L, 810L, 820L, 830L)); times.add(asTimes(1000, 0L, 0L, 0L)); times.add(asTimes(1000, 1L, 1L, 1L)); times.add(asTimes(1000, 999L, 999L, 999L)); times.add(asTimes(1000, 0L, 0L, 500L, 500L, 999L, 999L)); times.add(asTimes(1000, 250L, 250L, 500L, 500L, 750L, 750L)); times.add(asTimes(1000, 250L, 250L, 250L, 500L, 500L, 500L, 750L, 750L, 750L)); for (int i = 0; i < 10; i++) { times.add(generateTimes(rng, 10d)); } for (int i = 0; i < 10; i++) { times.add(generateTimes(rng, 30d)); } for (int i = 0; i < 5; i++) { final List<Double> ts = generateTimes(rng, 50d).list; final List<Double> newTs = newArrayList(); for (final double l : ts) { newTs.add(l); newTs.add(Math.min(999, Math.max(0, l + DoubleMath.roundToLong(rng.nextDouble() * 6d - 3d, RoundingMode.HALF_EVEN)))); } times.add(asTimesDouble(1000, newTs)); } for (int i = 0; i < 5; i++) { final List<Double> ts = generateTimes(rng, 100d).list; final List<Double> newTs = newArrayList(); System.out.println("num events: " + ts.size()); for (final double l : ts) { newTs.add(l); newTs.add(Math.min(999, Math.max(0, l + DoubleMath.roundToLong(rng.nextDouble() * 2d - 1d, RoundingMode.HALF_EVEN)))); newTs.add(Math.min(999, Math.max(0, l + DoubleMath.roundToLong(rng.nextDouble() * 2d - 1d, RoundingMode.HALF_EVEN)))); } times.add(asTimesDouble(1000, newTs)); } final List<Long> t2 = asList(10L, 30L, 50L, 70L, 90L); for (int i = 0; i < 5; i++) { final List<Long> c = newArrayList(); for (int j = 0; j < i + 1; j++) { c.addAll(t2); } Collections.sort(c); times.add(asTimes(100, c)); } final List<Long> t = asList(100L, 300L, 500L, 700L, 900L); for (int i = 0; i < 15; i++) { final List<Long> c = newArrayList(); for (int j = 0; j < i + 1; j++) { c.addAll(t); } Collections.sort(c); times.add(asTimes(1000, c)); } final List<Long> variant = newArrayList(); variant.addAll(t); for (int i = 0; i < 70; i++) { variant.add(100L); } Collections.sort(variant); times.add(asTimes(1000, variant)); checkState(variant.size() == 75); checkState(times.get(times.size() - 2).list.size() == 75, "", times.get(times.size() - 2).list.size()); for (int i = 0; i < 10; i++) { times.add(generateTimes(rng, (i + 1) * 100d)); } final ImmutableList<Long> all = ContiguousSet.create(Range.closedOpen(0L, 1000L), DiscreteDomain.longs()) .asList(); times.add(asTimes(1000, all)); final List<Long> more = newArrayList(all); for (final long l : all) { if (l % 2 == 0) { more.add(l); } } Collections.sort(more); times.add(asTimes(1000, more)); final List<Long> more2 = newArrayList(all); for (int i = 0; i < 200; i++) { more2.add(100L); } for (int i = 0; i < 100; i++) { more2.add(200L); } Collections.sort(more2); final List<Long> newMore = newArrayList(); for (int i = 0; i < more2.size(); i++) { newMore.add(more2.get(i) * 10L); } times.add(asTimes(1000, more2)); times.add(asTimes(10000, newMore)); for (int k = 0; k < 5; k++) { final List<Double> ts = newArrayList(generateTimes(rng, 20).list); final List<Double> additions = newArrayList(); for (int i = 0; i < ts.size(); i++) { if (i % 3 == 0) { for (int j = 0; j < k; j++) { additions.add(ts.get(i) + rng.nextDouble() * 10); } } } ts.addAll(additions); Collections.sort(ts); times.add(asTimesDouble(1000, ts)); } final Times regular = asTimes(10, 0L, 1L, 2L, 5L, 6L, 7L, 8L, 9L); for (int i = 1; i < 4; i++) { final List<Double> newList = newArrayList(); for (final double l : regular.list) { newList.add(l * Math.pow(10, i)); } times.add(asTimesDouble((int) (regular.length * Math.pow(10, i)), newList)); } times.add(asTimes(1000, 250L, 250L, 250L, 500L, 500L, 500L, 750L, 750L, 750L)); times.add(asTimes(1000, 250L, 500L, 500L, 500L, 500L, 500L, 500L, 500L, 750L)); times.add(asTimes(1000, 100L, 100L, 200L, 200L, 300L, 300L, 400L, 400L, 500L, 500L, 600L, 600L, 700L, 700L, 800L, 800L, 900L, 900L)); times.add(asTimes(1000, 100L, 200L, 200L, 200L, 200L, 200L, 200L, 200L, 200L, 200L, 200L, 300L, 400L, 500L, 600L, 700L, 800L, 900L)); times.add(asTimes(1000, 100L, 200L, 300L, 400L, 500L, 600L, 700L, 800L, 800L, 800L, 800L, 800L, 800L, 800L, 800L, 800L, 800L, 900L)); // times.subList(1, times.size()).clear(); for (int j = 0; j < ordersPerHour.length; j++) { System.out.println("=========" + ordersPerHour[j] + "========="); for (int i = 0; i < times.size(); i++) { System.out.println("----- " + i + " -----"); // System.out.println(times.get(i).length + " " + times.get(i).list); // final double dod2 = measureDynamismDistr(times.get(i).list, // times.get(i).length); final double dod8 = measureDynamism(times.get(i).list, times.get(i).length); // final double dod5 = measureDynamism2ndDerivative(times.get(i).list, // times.get(i).length); // final double dod6 = measureDynDeviationCount(times.get(i).list, // times.get(i).length); // final double dod7 = chi(times.get(i).list, // times.get(i).length); // System.out.println(dod); // System.out.printf("%1.3f%%\n", dod2 * 100d); System.out.printf("%1.3f%%\n", dod8 * 100d); // System.out.printf("%1.3f%%\n", dod5 * 100d); // System.out.printf("%1.3f%%\n", dod6 * 100d); // System.out.printf("%1.3f%%\n", dod7); final double name = Math.round(dod8 * 100000d) / 1000d; try { final File dest = new File("files/generator/times/orders" + Strings.padStart(Integer.toString(i), 3, '0') + "-" + name + ".times"); Files.createParentDirs(dest); Files.write(times.get(i).length + "\n" + Joiner.on("\n").join(times.get(i).list) + "\n", dest, Charsets.UTF_8); } catch (final IOException e) { throw new IllegalArgumentException(); } } } }