List of usage examples for org.apache.commons.math3.linear RealMatrix setColumn
void setColumn(int column, double[] array) throws OutOfRangeException, MatrixDimensionMismatchException;
From source file:hivemall.anomaly.SingularSpectrumTransform.java
@Override public void update(@Nonnull final Object arg, @Nonnull final double[] outScores) throws HiveException { double x = PrimitiveObjectInspectorUtils.getDouble(arg, oi); xRing.add(x).toArray(xSeries, true /* FIFO */); // need to wait until the buffer is filled if (!xRing.isFull()) { outScores[0] = 0.d;/*from w w w .j a va 2 s .c om*/ } else { // create past trajectory matrix and find its left singular vectors RealMatrix H = new Array2DRowRealMatrix(window, nPastWindow); for (int i = 0; i < nPastWindow; i++) { H.setColumn(i, Arrays.copyOfRange(xSeries, i, i + window)); } // create current trajectory matrix and find its left singular vectors RealMatrix G = new Array2DRowRealMatrix(window, nCurrentWindow); int currentHead = pastSize + currentOffset; for (int i = 0; i < nCurrentWindow; i++) { G.setColumn(i, Arrays.copyOfRange(xSeries, currentHead + i, currentHead + i + window)); } switch (scoreFunc) { case svd: outScores[0] = computeScoreSVD(H, G); break; case ika: outScores[0] = computeScoreIKA(H, G); break; default: throw new IllegalStateException("Unexpected score function: " + scoreFunc); } } }
From source file:com.caseystella.analytics.outlier.batch.rpca.AugmentedDickeyFuller.java
private void computeADFStatistics() { double[] y = diff(ts); RealMatrix designMatrix = null; int k = lag + 1; int n = ts.length - 1; RealMatrix z = MatrixUtils.createRealMatrix(laggedMatrix(y, k)); //has rows length(ts) - 1 - k + 1 RealVector zcol1 = z.getColumnVector(0); //has length length(ts) - 1 - k + 1 double[] xt1 = subsetArray(ts, k - 1, n - 1); //ts[k:(length(ts) - 1)], has length length(ts) - 1 - k + 1 double[] trend = sequence(k, n); //trend k:n, has length length(ts) - 1 - k + 1 if (k > 1) { RealMatrix yt1 = z.getSubMatrix(0, ts.length - 1 - k, 1, k - 1); //same as z but skips first column //build design matrix as cbind(xt1, 1, trend, yt1) designMatrix = MatrixUtils.createRealMatrix(ts.length - 1 - k + 1, 3 + k - 1); designMatrix.setColumn(0, xt1); designMatrix.setColumn(1, ones(ts.length - 1 - k + 1)); designMatrix.setColumn(2, trend); designMatrix.setSubMatrix(yt1.getData(), 0, 3); } else {/* ww w.j a v a 2s . c o m*/ //build design matrix as cbind(xt1, 1, tt) designMatrix = MatrixUtils.createRealMatrix(ts.length - 1 - k + 1, 3); designMatrix.setColumn(0, xt1); designMatrix.setColumn(1, ones(ts.length - 1 - k + 1)); designMatrix.setColumn(2, trend); } /*OLSMultipleLinearRegression regression = new OLSMultipleLinearRegression(); regression.setNoIntercept(true); regression.newSampleData(zcol1.toArray(), designMatrix.getData()); double[] beta = regression.estimateRegressionParameters(); double[] sd = regression.estimateRegressionParametersStandardErrors(); */ RidgeRegression regression = new RidgeRegression(designMatrix.getData(), zcol1.toArray()); regression.updateCoefficients(.0001); double[] beta = regression.getCoefficients(); double[] sd = regression.getStandarderrors(); double t = beta[0] / sd[0]; if (t <= PVALUE_THRESHOLD) { this.needsDiff = true; } else { this.needsDiff = false; } }
From source file:mase.app.allocation.AllocationProblem.java
private DescriptiveStatistics pairDistances(RealMatrix distMatrix) { DescriptiveStatistics ds = new DescriptiveStatistics(); distMatrix = distMatrix.copy();// ww w. j ava 2s . com for (int k = 0; k < numAgents; k++) { // find closest pair double min = Double.POSITIVE_INFINITY; int minI = -1, minJ = -1; for (int i = 0; i < numAgents; i++) { for (int j = 0; j < numAgents; j++) { double d = distMatrix.getEntry(i, j); if (!Double.isNaN(d) && d < min) { min = d; minI = i; minJ = j; } } } ds.addValue(min); distMatrix.setRow(minI, nulify); distMatrix.setColumn(minJ, nulify); } return ds; }
From source file:lirmm.inria.fr.math.BigSparseRealMatrixTest.java
@Test public void testSetColumn() { RealMatrix m = new BigSparseRealMatrix(subTestData); double[] mColumn3 = columnToArray(subColumn3); Assert.assertTrue(mColumn3[0] != m.getColumn(1)[0]); m.setColumn(1, mColumn3); checkArrays(mColumn3, m.getColumn(1)); try {// w ww .j a v a 2s .co m m.setColumn(-1, mColumn3); Assert.fail("Expecting OutOfRangeException"); } catch (OutOfRangeException ex) { // expected } try { m.setColumn(0, new double[5]); Assert.fail("Expecting MatrixDimensionMismatchException"); } catch (MatrixDimensionMismatchException ex) { // expected } }
From source file:org.orekit.frames.TransformTest.java
@Test public void testLinear() { RandomGenerator random = new Well19937a(0x14f6411217b148d8l); for (int n = 0; n < 100; ++n) { Transform t = randomTransform(random); // build an equivalent linear transform by extracting raw translation/rotation RealMatrix linearA = MatrixUtils.createRealMatrix(3, 4); linearA.setSubMatrix(t.getRotation().getMatrix(), 0, 0); Vector3D rt = t.getRotation().applyTo(t.getTranslation()); linearA.setEntry(0, 3, rt.getX()); linearA.setEntry(1, 3, rt.getY()); linearA.setEntry(2, 3, rt.getZ()); // build an equivalent linear transform by observing transformed points RealMatrix linearB = MatrixUtils.createRealMatrix(3, 4); Vector3D p0 = t.transformPosition(Vector3D.ZERO); Vector3D pI = t.transformPosition(Vector3D.PLUS_I).subtract(p0); Vector3D pJ = t.transformPosition(Vector3D.PLUS_J).subtract(p0); Vector3D pK = t.transformPosition(Vector3D.PLUS_K).subtract(p0); linearB.setColumn(0, new double[] { pI.getX(), pI.getY(), pI.getZ() }); linearB.setColumn(1, new double[] { pJ.getX(), pJ.getY(), pJ.getZ() }); linearB.setColumn(2, new double[] { pK.getX(), pK.getY(), pK.getZ() }); linearB.setColumn(3, new double[] { p0.getX(), p0.getY(), p0.getZ() }); // both linear transforms should be equal Assert.assertEquals(0.0, linearB.subtract(linearA).getNorm(), 1.0e-15 * linearA.getNorm()); for (int i = 0; i < 100; ++i) { Vector3D p = randomVector(1.0e3, random); Vector3D q = t.transformPosition(p); double[] qA = linearA.operate(new double[] { p.getX(), p.getY(), p.getZ(), 1.0 }); Assert.assertEquals(q.getX(), qA[0], 1.0e-13 * p.getNorm()); Assert.assertEquals(q.getY(), qA[1], 1.0e-13 * p.getNorm()); Assert.assertEquals(q.getZ(), qA[2], 1.0e-13 * p.getNorm()); double[] qB = linearB.operate(new double[] { p.getX(), p.getY(), p.getZ(), 1.0 }); Assert.assertEquals(q.getX(), qB[0], 1.0e-10 * p.getNorm()); Assert.assertEquals(q.getY(), qB[1], 1.0e-10 * p.getNorm()); Assert.assertEquals(q.getZ(), qB[2], 1.0e-10 * p.getNorm()); }//from w w w . j a v a 2s .c o m } }
From source file:put.ci.cevo.framework.algorithms.ApacheCMAES.java
/** * {@inheritDoc}/*w ww.j a v a 2 s . co m*/ */ @Override protected PointValuePair doOptimize() { // -------------------- Initialization -------------------------------- isMinimize = getGoalType().equals(GoalType.MINIMIZE); final double[] guess = getStartPoint(); // number of objective variables/problem dimension dimension = guess.length; initializeCMA(guess); iterations = 0; double bestValue = (isMinimize ? Double.MAX_VALUE : Double.MIN_VALUE); push(fitnessHistory, bestValue); PointValuePair optimum = new PointValuePair(getStartPoint(), isMinimize ? bestValue : -bestValue); PointValuePair lastResult = null; // -------------------- Generation Loop -------------------------------- EvaluatedPopulation<double[]> evaluatedPopulation = null; Stopwatch stopwatch = Stopwatch.createUnstarted(); generationLoop: for (iterations = 1; iterations <= maxIterations; iterations++) { stopwatch.reset(); stopwatch.start(); incrementIterationCount(); // Generate and evaluate lambda offspring final RealMatrix arz = randn1(dimension, lambda); final RealMatrix arx = zeros(dimension, lambda); final double[] fitness = new double[lambda]; // generate random offspring for (int k = 0; k < lambda; k++) { RealMatrix arxk = null; for (int i = 0; i < checkFeasableCount + 1; i++) { if (diagonalOnly <= 0) { arxk = xmean.add(BD.multiply(arz.getColumnMatrix(k)).scalarMultiply(sigma)); // m + sig * Normal(0,C) } else { arxk = xmean.add(times(diagD, arz.getColumnMatrix(k)).scalarMultiply(sigma)); } //if (i >= checkFeasableCount || // fitfun.isFeasible(arxk.getColumn(0))) { // break; //} // regenerate random arguments for row arz.setColumn(k, randn(dimension)); } copyColumn(arxk, 0, arx, k); //try { // valuePenaltyPairs[k] = fitfun.value(arx.getColumn(k)); // compute fitness //} catch (TooManyEvaluationsException e) { // break generationLoop; //} } double newPopTime = stopwatch.elapsed(TimeUnit.MILLISECONDS) / 1000.0; stopwatch.reset(); stopwatch.start(); ArrayList<double[]> population = new ArrayList<>(lambda); // This is mine. I ignore constraints. for (int k = 0; k < lambda; ++k) { population.add(arx.getColumn(k)); } evaluatedPopulation = populationEvaluator.evaluate(population, iterations - 1, random); final ValuePenaltyPair[] valuePenaltyPairs = new ValuePenaltyPair[lambda]; for (int k = 0; k < lambda; ++k) { valuePenaltyPairs[k] = new ValuePenaltyPair(evaluatedPopulation.getPopulation().get(k).getFitness(), 0.0); } // Compute fitnesses by adding value and penalty after scaling by value range. double valueRange = valueRange(valuePenaltyPairs); for (int iValue = 0; iValue < valuePenaltyPairs.length; iValue++) { fitness[iValue] = valuePenaltyPairs[iValue].value + valuePenaltyPairs[iValue].penalty * valueRange; if (!isMinimize) fitness[iValue] = -fitness[iValue]; } double evalTime = stopwatch.elapsed(TimeUnit.MILLISECONDS) / 1000.0; stopwatch.reset(); stopwatch.start(); // Sort by fitness and compute weighted mean into xmean final int[] arindex = sortedIndices(fitness); // Calculate new xmean, this is selection and recombination final RealMatrix xold = xmean; // for speed up of Eq. (2) and (3) final RealMatrix bestArx = selectColumns(arx, MathArrays.copyOf(arindex, mu)); xmean = bestArx.multiply(weights); final RealMatrix bestArz = selectColumns(arz, MathArrays.copyOf(arindex, mu)); final RealMatrix zmean = bestArz.multiply(weights); final boolean hsig = updateEvolutionPaths(zmean, xold); if (diagonalOnly <= 0) { updateCovariance(hsig, bestArx, arz, arindex, xold); } else { updateCovarianceDiagonalOnly(hsig, bestArz); } // Adapt step size sigma - Eq. (5) sigma *= FastMath.exp(FastMath.min(1, (normps / chiN - 1) * cs / damps)); final double bestFitness = fitness[arindex[0]]; final double worstFitness = fitness[arindex[arindex.length - 1]]; if (bestValue > bestFitness) { bestValue = bestFitness; lastResult = optimum; optimum = new PointValuePair(bestArx.getColumn(0), isMinimize ? bestFitness : -bestFitness); if (getConvergenceChecker() != null && lastResult != null && getConvergenceChecker().converged(iterations, optimum, lastResult)) { break generationLoop; } } // handle termination criteria // Break, if fitness is good enough if (stopFitness != 0 && bestFitness < (isMinimize ? stopFitness : -stopFitness)) { break generationLoop; } final double[] sqrtDiagC = sqrt(diagC).getColumn(0); final double[] pcCol = pc.getColumn(0); for (int i = 0; i < dimension; i++) { if (sigma * FastMath.max(FastMath.abs(pcCol[i]), sqrtDiagC[i]) > stopTolX) { break; } if (i >= dimension - 1) { break generationLoop; } } for (int i = 0; i < dimension; i++) { if (sigma * sqrtDiagC[i] > stopTolUpX) { break generationLoop; } } final double historyBest = min(fitnessHistory); final double historyWorst = max(fitnessHistory); if (iterations > 2 && FastMath.max(historyWorst, worstFitness) - FastMath.min(historyBest, bestFitness) < stopTolFun) { break generationLoop; } if (iterations > fitnessHistory.length && historyWorst - historyBest < stopTolHistFun) { break generationLoop; } // condition number of the covariance matrix exceeds 1e14 if (max(diagD) / min(diagD) > 1e7) { break generationLoop; } // user defined termination if (getConvergenceChecker() != null) { final PointValuePair current = new PointValuePair(bestArx.getColumn(0), isMinimize ? bestFitness : -bestFitness); if (lastResult != null && getConvergenceChecker().converged(iterations, current, lastResult)) { break generationLoop; } lastResult = current; } // Adjust step size in case of equal function values (flat fitness) if (bestValue == fitness[arindex[(int) (0.1 + lambda / 4.)]]) { sigma *= FastMath.exp(0.2 + cs / damps); } if (iterations > 2 && FastMath.max(historyWorst, bestFitness) - FastMath.min(historyBest, bestFitness) == 0) { sigma *= FastMath.exp(0.2 + cs / damps); } // store best in history push(fitnessHistory, bestFitness); if (generateStatistics) { statisticsSigmaHistory.add(sigma); statisticsFitnessHistory.add(bestFitness); statisticsMeanHistory.add(xmean.transpose()); statisticsDHistory.add(diagD.transpose().scalarMultiply(1E5)); } double cmaesTime = stopwatch.elapsed(TimeUnit.MILLISECONDS) / 1000.0; stopwatch.reset(); stopwatch.start(); listener.onNextIteraction(evaluatedPopulation); double listernerTime = stopwatch.elapsed(TimeUnit.MILLISECONDS) / 1000.0; logger.info(String.format("NewPop: %.2f, Eval: %.2f, CMAES: %.2f, Listerner: %.2f", newPopTime, evalTime, cmaesTime, listernerTime)); } listener.onLastIteraction(evaluatedPopulation); return optimum; }
From source file:stats.SpearmansCorrelation.java
/** * Applies rank transform to each of the columns of <code>matrix</code> using * the current <code>rankingAlgorithm</code>. * * @param matrix//from ww w.j a va2s . c om * matrix to transform * @return a rank-transformed matrix */ private RealMatrix rankTransform(final RealMatrix matrix) { RealMatrix transformed = null; if (rankingAlgorithm instanceof NaturalRanking && ((NaturalRanking) rankingAlgorithm).getNanStrategy() == NaNStrategy.REMOVED) { final Set<Integer> nanPositions = new HashSet<Integer>(); for (int i = 0; i < matrix.getColumnDimension(); i++) { nanPositions.addAll(getNaNPositions(matrix.getColumn(i))); } // if we have found NaN values, we have to update the matrix size if (!nanPositions.isEmpty()) { transformed = new BlockRealMatrix(matrix.getRowDimension() - nanPositions.size(), matrix.getColumnDimension()); for (int i = 0; i < transformed.getColumnDimension(); i++) { transformed.setColumn(i, removeValues(matrix.getColumn(i), nanPositions)); } } } if (transformed == null) { transformed = matrix.copy(); } for (int i = 0; i < transformed.getColumnDimension(); i++) { transformed.setColumn(i, rankingAlgorithm.rank(transformed.getColumn(i))); } return transformed; }