List of usage examples for org.apache.commons.math3.linear RealMatrix scalarMultiply
RealMatrix scalarMultiply(double d);
From source file:org.rhwlab.variationalbayesian.GaussianMixture.java
public void Mstep() { alphaBetaNu();/*from www . ja v a2 s. c o m*/ // compute the means of the components for (int k = 0; k < K; ++k) { if (inModel[k]) { RealVector v1 = xBar[k].mapMultiply(N[k]); RealVector v2 = mu0.mapMultiply(beta0); RealVector v3 = v2.add(v1); m[k] = v3.mapDivide(beta[k]); } } // compute the precision matrices for (int k = 0; k < K; ++k) { if (inModel[k]) { RealVector del = xBar[k].subtract(mu0); RealMatrix del2 = del.outerProduct(del); double f = beta0 * N[k] / (beta0 + N[k]); RealMatrix mat = del2.scalarMultiply(f); RealMatrix NS = S[k].scalarMultiply(N[k]); Winverse[k] = W0inverse.add(NS).add(mat); LUDecomposition cd = new LUDecomposition(Winverse[k]); W[k] = cd.getSolver().getInverse(); detWinv[k] = cd.getDeterminant(); detW[k] = 1.0 / cd.getDeterminant(); } } lambdaTilde(); }
From source file:org.rhwlab.variationalbayesian.OldFaithfulDataSource.java
static public void main(String[] args) throws Exception { OldFaithfulDataSource source = new OldFaithfulDataSource( "/net/waterston/vol2/home/gevirl/Downloads/VBEMGMM/faithful.txt"); GaussianMixture gm = new GaussianMixture(); gm.setSource(source);/*w w w.ja v a2 s. c o m*/ gm.setAlpha0(0.001); gm.setBeta0(1.0); gm.setNu0(20.0); RealMatrix W0 = MatrixUtils.createRealIdentityMatrix(source.getD()); W0 = W0.scalarMultiply(200.0); gm.setW0(W0); gm.init(15); gm.run(); }
From source file:org.rhwlab.variationalbayesian.SuperVoxelGaussianMixture.java
public void Mstep() { lnPi();/*from www. j av a 2 s .c o m*/ alphaBetaNuLambda(); // compute the means of the components for (int k = 0; k < K; ++k) { RealVector v1 = xBar[k].mapMultiply(N[k]); RealVector v2 = mu0.mapMultiply(beta0); RealVector v3 = v2.add(v1); m[k] = v3.mapDivide(beta[k]); } // compute the precision matrices for (int k = 0; k < K; ++k) { RealVector del = xBar[k].subtract(mu0); RealMatrix del2 = del.outerProduct(del); double f = beta0 * N[k] / (beta0 + N[k]); RealMatrix mat = del2.scalarMultiply(f); RealMatrix NS = S[k].scalarMultiply(N[k]); RealMatrix Winverse = W0inverse.add(NS).add(mat); LUDecomposition cd = new LUDecomposition(Winverse); W[k] = cd.getSolver().getInverse(); detW[k] = 1.0 / cd.getDeterminant(); int auiosdfu = 0; } }
From source file:org.um.feri.ears.algorithms.moo.dbea.DBEA.java
/** * Updates the ideal point and intercepts given the new solution. * /* w w w .j a v a 2 s . c om*/ * @param solution the new solution */ void updateIdealPointAndIntercepts(MOSolutionBase<Type> solution) { if (!solution.violatesConstraints()) { // update the ideal point for (int j = 0; j < num_obj; j++) { idealPoint[j] = Math.min(idealPoint[j], solution.getObjective(j)); intercepts[j] = Math.max(intercepts[j], solution.getObjective(j)); } // compute the axis intercepts ParetoSolution<Type> feasibleSolutions = getFeasibleSolutions(population); feasibleSolutions.add(solution); ParetoSolution<Type> nondominatedSolutions = getNondominatedFront(feasibleSolutions); if (!nondominatedSolutions.isEmpty()) { // find the points with the largest value in each objective ParetoSolution<Type> extremePoints = new ParetoSolution<Type>(); for (int i = 0; i < num_obj; i++) { extremePoints.add(largestObjectiveValue(i, nondominatedSolutions)); } if (numberOfUniqueSolutions(extremePoints) != num_obj) { for (int i = 0; i < num_obj; i++) { intercepts[i] = extremePoints.get(i).getObjective(i); } } else { try { RealMatrix b = new Array2DRowRealMatrix(num_obj, 1); RealMatrix A = new Array2DRowRealMatrix(num_obj, num_obj); for (int i = 0; i < num_obj; i++) { b.setEntry(i, 0, 1.0); for (int j = 0; j < num_obj; j++) { A.setEntry(i, j, extremePoints.get(i).getObjective(j)); } } double numerator = new LUDecomposition(A).getDeterminant(); b.scalarMultiply(numerator); RealMatrix normal = MatrixUtils.inverse(A).multiply(b); for (int i = 0; i < num_obj; i++) { intercepts[i] = numerator / normal.getEntry(i, 0); if (intercepts[i] <= 0 || Double.isNaN(intercepts[i]) || Double.isInfinite(intercepts[i])) { intercepts[i] = extremePoints.get(i).getObjective(i); } } } catch (RuntimeException e) { for (int i = 0; i < num_obj; i++) { intercepts[i] = extremePoints.get(i).getObjective(i); } } } } } }
From source file:pl.matrix.core.MatrixCalculator.java
public Matrix multiplyByScalar(Matrix a, double s) { RealMatrix aRealMatrix = a.toRealMatrix(); Matrix result = new Matrix(aRealMatrix.scalarMultiply(s)); return result; }
From source file:playground.sergioo.facilitiesGenerator2012.WorkFacilitiesGeneration.java
private static Set<PointPerson> getPCATransformation(Collection<PointPerson> points) { RealMatrix pointsM = new Array2DRowRealMatrix(points.iterator().next().getDimension(), points.size()); int k = 0;/*from w w w.j a v a2 s. com*/ for (PointND<Double> point : points) { for (int f = 0; f < point.getDimension(); f++) pointsM.setEntry(f, k, point.getElement(f)); k++; } RealMatrix means = new Array2DRowRealMatrix(pointsM.getRowDimension(), 1); for (int r = 0; r < means.getRowDimension(); r++) { double mean = 0; for (int c = 0; c < pointsM.getColumnDimension(); c++) mean += pointsM.getEntry(r, c) / pointsM.getColumnDimension(); means.setEntry(r, 0, mean); } RealMatrix deviations = new Array2DRowRealMatrix(pointsM.getRowDimension(), pointsM.getColumnDimension()); for (int r = 0; r < deviations.getRowDimension(); r++) for (int c = 0; c < deviations.getColumnDimension(); c++) deviations.setEntry(r, c, pointsM.getEntry(r, c) - means.getEntry(r, 0)); RealMatrix covariance = deviations.multiply(deviations.transpose()) .scalarMultiply(1 / (double) pointsM.getColumnDimension()); EigenDecomposition eigenDecomposition = new EigenDecomposition(covariance, 0); RealMatrix eigenVectorsT = eigenDecomposition.getVT(); RealVector eigenValues = new ArrayRealVector(eigenDecomposition.getD().getRowDimension()); for (int r = 0; r < eigenDecomposition.getD().getRowDimension(); r++) eigenValues.setEntry(r, eigenDecomposition.getD().getEntry(r, r)); for (int i = 0; i < eigenValues.getDimension(); i++) { for (int j = i + 1; j < eigenValues.getDimension(); j++) if (eigenValues.getEntry(i) < eigenValues.getEntry(j)) { double tempValue = eigenValues.getEntry(i); eigenValues.setEntry(i, eigenValues.getEntry(j)); eigenValues.setEntry(j, tempValue); RealVector tempVector = eigenVectorsT.getRowVector(i); eigenVectorsT.setRowVector(i, eigenVectorsT.getRowVector(j)); eigenVectorsT.setRowVector(j, tempVector); } eigenVectorsT.setRowVector(i, eigenVectorsT.getRowVector(i).mapMultiply(Math.sqrt(1 / eigenValues.getEntry(i)))); } RealVector standardDeviations = new ArrayRealVector(pointsM.getRowDimension()); for (int r = 0; r < covariance.getRowDimension(); r++) standardDeviations.setEntry(r, Math.sqrt(covariance.getEntry(r, r))); double zValue = standardDeviations.dotProduct(new ArrayRealVector(pointsM.getRowDimension(), 1)); RealMatrix zScore = deviations.scalarMultiply(1 / zValue); pointsM = eigenVectorsT.multiply(zScore); Set<PointPerson> pointsC = new HashSet<PointPerson>(); k = 0; for (PointPerson point : points) { PointPerson pointC = new PointPerson(point.getId(), point.getOccupation(), new Double[] { pointsM.getEntry(0, k), pointsM.getEntry(1, k) }, point.getPlaceType()); pointC.setWeight(point.getWeight()); pointsC.add(pointC); k++; } return pointsC; }
From source file:put.ci.cevo.framework.algorithms.ApacheCMAES.java
/** * Initialization of the dynamic search parameters * * @param guess Initial guess for the arguments of the fitness function. *///from ww w. j a va2s. com private void initializeCMA(double[] guess) { if (lambda <= 0) { throw new NotStrictlyPositiveException(lambda); } // initialize sigma final double[][] sigmaArray = new double[guess.length][1]; for (int i = 0; i < guess.length; i++) { sigmaArray[i][0] = inputSigma[i]; } final RealMatrix insigma = new Array2DRowRealMatrix(sigmaArray, false); sigma = max(insigma); // overall standard deviation // initialize termination criteria stopTolUpX = 1e3 * max(insigma); stopTolX = 1e-11 * max(insigma); stopTolFun = 1e-12; stopTolHistFun = 1e-13; // initialize selection strategy parameters mu = lambda / 2; // number of parents/points for recombination logMu2 = FastMath.log(mu + 0.5); weights = log(sequence(1, mu, 1)).scalarMultiply(-1).scalarAdd(logMu2); double sumw = 0; double sumwq = 0; for (int i = 0; i < mu; i++) { double w = weights.getEntry(i, 0); sumw += w; sumwq += w * w; } weights = weights.scalarMultiply(1 / sumw); mueff = sumw * sumw / sumwq; // variance-effectiveness of sum w_i x_i // initialize dynamic strategy parameters and constants cc = (4 + mueff / dimension) / (dimension + 4 + 2 * mueff / dimension); cs = (mueff + 2) / (dimension + mueff + 3.); damps = (1 + 2 * FastMath.max(0, FastMath.sqrt((mueff - 1) / (dimension + 1)) - 1)) * FastMath.max(0.3, 1 - dimension / (1e-6 + maxIterations)) + cs; // minor increment ccov1 = 2 / ((dimension + 1.3) * (dimension + 1.3) + mueff); ccovmu = FastMath.min(1 - ccov1, 2 * (mueff - 2 + 1 / mueff) / ((dimension + 2) * (dimension + 2) + mueff)); ccov1Sep = FastMath.min(1, ccov1 * (dimension + 1.5) / 3); ccovmuSep = FastMath.min(1 - ccov1, ccovmu * (dimension + 1.5) / 3); chiN = FastMath.sqrt(dimension) * (1 - 1 / ((double) 4 * dimension) + 1 / ((double) 21 * dimension * dimension)); // intialize CMA internal values - updated each generation xmean = MatrixUtils.createColumnRealMatrix(guess); // objective variables diagD = insigma.scalarMultiply(1 / sigma); diagC = square(diagD); pc = zeros(dimension, 1); // evolution paths for C and sigma ps = zeros(dimension, 1); // B defines the coordinate system normps = ps.getFrobeniusNorm(); B = eye(dimension, dimension); D = ones(dimension, 1); // diagonal D defines the scaling BD = times(B, repmat(diagD.transpose(), dimension, 1)); C = eye(dimension, dimension); /// WJ: originally it was: //C = B.multiply(diag(square(D)).multiply(B.transpose())); // covariance historySize = 10 + (int) (3 * 10 * dimension / (double) lambda); fitnessHistory = new double[historySize]; // history of fitness values for (int i = 0; i < historySize; i++) { fitnessHistory[i] = Double.MAX_VALUE; } }
From source file:put.ci.cevo.framework.algorithms.ApacheCMAES.java
/** * Update of the covariance matrix C./* w ww . j a v a 2 s . c o m*/ * * @param hsig Flag indicating a small correction. * @param bestArx Fitness-sorted matrix of the argument vectors producing the current offspring. * @param arz Unsorted matrix containing the gaussian random values of the current offspring. * @param arindex Indices indicating the fitness-order of the current offspring. * @param xold xmean matrix of the previous generation. */ private void updateCovariance(boolean hsig, final RealMatrix bestArx, final RealMatrix arz, final int[] arindex, final RealMatrix xold) { double negccov = 0; if (ccov1 + ccovmu > 0) { final RealMatrix arpos = bestArx.subtract(repmat(xold, 1, mu)).scalarMultiply(1 / sigma); // mu difference vectors final RealMatrix roneu = pc.multiply(pc.transpose()).scalarMultiply(ccov1); // rank one update // minor correction if hsig==false double oldFac = hsig ? 0 : ccov1 * cc * (2 - cc); oldFac += 1 - ccov1 - ccovmu; if (isActiveCMA) { // Adapt covariance matrix C active CMA negccov = (1 - ccovmu) * 0.25 * mueff / (FastMath.pow(dimension + 2, 1.5) + 2 * mueff); // keep at least 0.66 in all directions, small popsize are most // critical final double negminresidualvariance = 0.66; // where to make up for the variance loss final double negalphaold = 0.5; // prepare vectors, compute negative updating matrix Cneg final int[] arReverseIndex = reverse(arindex); RealMatrix arzneg = selectColumns(arz, MathArrays.copyOf(arReverseIndex, mu)); RealMatrix arnorms = sqrt(sumRows(square(arzneg))); final int[] idxnorms = sortedIndices(arnorms.getRow(0)); final RealMatrix arnormsSorted = selectColumns(arnorms, idxnorms); final int[] idxReverse = reverse(idxnorms); final RealMatrix arnormsReverse = selectColumns(arnorms, idxReverse); arnorms = divide(arnormsReverse, arnormsSorted); final int[] idxInv = inverse(idxnorms); final RealMatrix arnormsInv = selectColumns(arnorms, idxInv); // check and set learning rate negccov final double negcovMax = (1 - negminresidualvariance) / square(arnormsInv).multiply(weights).getEntry(0, 0); if (negccov > negcovMax) { negccov = negcovMax; } arzneg = times(arzneg, repmat(arnormsInv, dimension, 1)); final RealMatrix artmp = BD.multiply(arzneg); final RealMatrix Cneg = artmp.multiply(diag(weights)).multiply(artmp.transpose()); oldFac += negalphaold * negccov; C = C.scalarMultiply(oldFac).add(roneu) // regard old matrix .add(arpos.scalarMultiply( // plus rank one update ccovmu + (1 - negalphaold) * negccov) // plus rank mu update .multiply(times(repmat(weights, 1, dimension), arpos.transpose()))) .subtract(Cneg.scalarMultiply(negccov)); } else { // Adapt covariance matrix C - nonactive C = C.scalarMultiply(oldFac) // regard old matrix .add(roneu) // plus rank one update .add(arpos.scalarMultiply(ccovmu) // plus rank mu update .multiply(times(repmat(weights, 1, dimension), arpos.transpose()))); } } updateBDFast(negccov); }