List of usage examples for org.apache.commons.math3.linear RealVector mapMultiply
public RealVector mapMultiply(double d)
From source file:edu.stanford.cfuller.imageanalysistools.fitting.BisquareLinearFit.java
/** * Calculates the leverages of data points for least squares fitting (assuming equal variances). * * @param indVarValues The values of the independent variable used for the fitting. * @return a RealVector containing a leverage value for each independent variable value. *//*from w ww . j a va 2 s . c om*/ protected RealVector calculateLeverages(RealVector indVarValues) { RealMatrix indVarMatrix = null; if (this.noIntercept) { indVarMatrix = new Array2DRowRealMatrix(indVarValues.getDimension(), 1); } else { indVarMatrix = new Array2DRowRealMatrix(indVarValues.getDimension(), 2); } indVarMatrix.setColumnVector(0, indVarValues); if (!this.noIntercept) { indVarMatrix.setColumnVector(1, indVarValues.mapMultiply(0).mapAdd(1)); } RealVector leverages = new ArrayRealVector(indVarValues.getDimension()); QRDecomposition xQR = new QRDecomposition(indVarMatrix); RealMatrix xR = xQR.getR(); int smallerDim = xR.getRowDimension() < xR.getColumnDimension() ? xR.getRowDimension() : xR.getColumnDimension(); RealMatrix xRSq = xR.getSubMatrix(0, smallerDim - 1, 0, smallerDim - 1); QRDecomposition xRQR = new QRDecomposition(xRSq); RealMatrix xRInv = xRQR.getSolver().getInverse(); RealMatrix xxRInv = indVarMatrix.multiply(xRInv); for (int i = 0; i < indVarValues.getDimension(); i++) { double sum = 0; for (int j = 0; j < xxRInv.getColumnDimension(); j++) { sum += Math.pow(xxRInv.getEntry(i, j), 2); } leverages.setEntry(i, sum); } return leverages; }
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); }/*www . ja va 2 s . 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:edu.utexas.cs.tactex.tariffoptimization.OptimizerWrapperGradientAscent.java
@Override public TreeMap<Double, TariffSpecification> findOptimum(TariffUtilityEstimate tariffUtilityEstimate, int NUM_RATES, int numEval) { double[] startingVertex = new double[NUM_RATES]; // start from the fixed-rate tariff's offset Arrays.fill(startingVertex, 0.0); // temporary solution - getting fixed rate tariff to determine step size and scaling STEP_SIZE proportionally TariffSpecification bestFixedRateSpec = tariffUtilityEstimate.getCorrespondingSpec(startingVertex); double bestFixedRate = bestFixedRateSpec.getRates().get(0).getValue(); double rateRatio = bestFixedRate / REFERENCE_RATE; // Note: using rateRatio has not been good enough, so trying powers of it (remember it's > 1) //STEP_SIZE = Math.max(REFERENCE_STEP_SIZE, (rateRatio * rateRatio) * REFERENCE_STEP_SIZE); STEP_SIZE = Math.max(REFERENCE_STEP_SIZE, Math.abs(rateRatio) * REFERENCE_STEP_SIZE); log.debug("STEP_SIZE = " + STEP_SIZE + " REFERENCE_STEP_SIZE=" + REFERENCE_STEP_SIZE + " bestFixedRate=" + bestFixedRate + " REFERENCE_RATE=" + REFERENCE_RATE); evaluations = 0;/*from www. j av a2s . co m*/ TreeMap<Double, TariffSpecification> eval2TOUTariff = new TreeMap<Double, TariffSpecification>(); // OUTER LOOP //for( STEP_SIZE = 0.005; STEP_SIZE < 0.100; STEP_SIZE += 0.005) { // log.info("STARTING A LOOP: STEP_SIZE=" + STEP_SIZE); // first compute numerical gradient RealVector gradient = new ArrayRealVector(NUM_RATES); for (int i = 0; i < NUM_RATES; ++i) { gradient.setEntry(i, computePartialDerivative(i, tariffUtilityEstimate, NUM_RATES, eval2TOUTariff)); } gradient = gradient.unitVector(); // taking steps in the gradient direction double previousPointValue = -Double.MAX_VALUE; final double alpha = STEP_SIZE; RealVector rateOffsets = new ArrayRealVector(NUM_RATES); // initializes with 0s? double currentPointValue = evaluatePoint(tariffUtilityEstimate, rateOffsets.toArray()); eval2TOUTariff.put(currentPointValue, tariffUtilityEstimate.getCorrespondingSpec(rateOffsets.toArray())); while (!converged(currentPointValue, previousPointValue) && evaluations < MAX_EVALUATIONS) { previousPointValue = currentPointValue; rateOffsets = rateOffsets.add(gradient.mapMultiply(alpha)); currentPointValue = evaluatePoint(tariffUtilityEstimate, rateOffsets.toArray()); eval2TOUTariff.put(currentPointValue, tariffUtilityEstimate.getCorrespondingSpec(rateOffsets.toArray())); } log.info("gradient ascent finished after " + evaluations + " evaluations"); //} // return map return eval2TOUTariff; }
From source file:com.joptimizer.optimizers.JOptimizerTest.java
/** * Quadratic objective, no constraints./*from www .j a va2s. c o m*/ */ public void testNewtownUnconstrained() throws Exception { Log.d(MainActivity.JOPTIMIZER_LOGTAG, "testNewtownUnconstrained"); RealMatrix PMatrix = new Array2DRowRealMatrix( new double[][] { { 1.68, 0.34, 0.38 }, { 0.34, 3.09, -1.59 }, { 0.38, -1.59, 1.54 } }); RealVector qVector = new ArrayRealVector(new double[] { 0.018, 0.025, 0.01 }); // Objective function. double theta = 0.01522; RealMatrix P = PMatrix.scalarMultiply(theta); RealVector q = qVector.mapMultiply(-1); PDQuadraticMultivariateRealFunction objectiveFunction = new PDQuadraticMultivariateRealFunction(P.getData(), q.toArray(), 0); OptimizationRequest or = new OptimizationRequest(); or.setF0(objectiveFunction); or.setInitialPoint(new double[] { 0.04, 0.50, 0.46 }); //optimization JOptimizer opt = new JOptimizer(); opt.setOptimizationRequest(or); int returnCode = opt.optimize(); if (returnCode == OptimizationResponse.FAILED) { fail(); } OptimizationResponse response = opt.getOptimizationResponse(); double[] sol = response.getSolution(); Log.d(MainActivity.JOPTIMIZER_LOGTAG, "sol : " + ArrayUtils.toString(sol)); Log.d(MainActivity.JOPTIMIZER_LOGTAG, "value : " + objectiveFunction.value(sol)); // we already know the analytic solution of the problem // sol = -invQ * C RealMatrix QInv = Utils.squareMatrixInverse(P); RealVector benchSol = QInv.operate(q).mapMultiply(-1); Log.d(MainActivity.JOPTIMIZER_LOGTAG, "benchSol : " + ArrayUtils.toString(benchSol.toArray())); Log.d(MainActivity.JOPTIMIZER_LOGTAG, "benchValue : " + objectiveFunction.value(benchSol.toArray())); assertEquals(benchSol.getEntry(0), sol[0], 0.000000000000001); assertEquals(benchSol.getEntry(1), sol[1], 0.000000000000001); assertEquals(benchSol.getEntry(2), sol[2], 0.000000000000001); }
From source file:edu.stanford.cfuller.colocalization3d.Colocalization3DMain.java
private String formatPositionData(java.util.List<ImageObject> imageObjects, Correction c) { StringBuilder sb = new StringBuilder(); StringBuilder sbSingle = new StringBuilder(); String currentImageName = ""; File lastFile = null;/*from ww w . ja v a2s. c om*/ for (int o = 0; o < imageObjects.size(); o++) { ImageObject currObject = imageObjects.get(o); String objectImageName = currObject.getImageID(); if (!currentImageName.equals(objectImageName)) { if (!(sbSingle.length() == 0)) { sb.append(sbSingle.toString()); sbSingle = new StringBuilder(); } sb.append(objectImageName + "\n"); currentImageName = objectImageName; } RealVector correction = null; try { correction = c.correctPosition( currObject.getPositionForChannel(c.getReferenceChannelIndex()).getEntry(0), currObject.getPositionForChannel(c.getReferenceChannelIndex()).getEntry(1)); } catch (UnableToCorrectException e) { continue; } sbSingle.append(currObject.getLabel()); sbSingle.append(" "); sbSingle.append(currObject.getPositionForChannel(c.getReferenceChannelIndex()).getEntry(0)); sbSingle.append(" "); sbSingle.append(currObject.getPositionForChannel(c.getReferenceChannelIndex()).getEntry(1)); sbSingle.append(" "); sbSingle.append(currObject.getPositionForChannel(c.getReferenceChannelIndex()).getEntry(2)); sbSingle.append(" "); sbSingle.append(currObject.getPositionForChannel(c.getCorrectionChannelIndex()).getEntry(0)); sbSingle.append(" "); sbSingle.append(currObject.getPositionForChannel(c.getCorrectionChannelIndex()).getEntry(1)); sbSingle.append(" "); sbSingle.append(currObject.getPositionForChannel(c.getCorrectionChannelIndex()).getEntry(2)); sbSingle.append(" "); RealVector corrPos = null; RealVector modCorrection = correction; if (this.parameters.hasKeyAndTrue("flip_channels_at_end")) { modCorrection = correction.mapMultiply(-1.0); } if (this.parameters.hasKeyAndTrue("inverted_z_axis")) { modCorrection = correction.mapMultiply(1.0); modCorrection.setEntry(2, -1.0 * modCorrection.getEntry(2)); } corrPos = currObject.getPositionForChannel(c.getCorrectionChannelIndex()).subtract(modCorrection); sbSingle.append(corrPos.getEntry(0)); sbSingle.append(" "); sbSingle.append(corrPos.getEntry(1)); sbSingle.append(" "); sbSingle.append(corrPos.getEntry(2)); sbSingle.append(" "); sbSingle.append( currObject.getFitParametersByChannel().get(c.getReferenceChannelIndex()).getAmplitude()); sbSingle.append(" "); sbSingle.append( currObject.getFitParametersByChannel().get(c.getCorrectionChannelIndex()).getAmplitude()); sbSingle.append("\n"); } sb.append(sbSingle.toString()); return sb.toString(); }
From source file:org.grouplens.samantha.modeler.svdfeature.SVDFeature.java
public List<StochasticOracle> getStochasticOracle(List<LearningInstance> instances) { List<StochasticOracle> oracles = new ArrayList<>(instances.size()); for (LearningInstance inIns : instances) { SVDFeatureInstance ins = (SVDFeatureInstance) inIns; StochasticOracle orc = new StochasticOracle(); RealVector ufactSum = MatrixUtils.createRealVector(new double[factDim]); RealVector ifactSum = MatrixUtils.createRealVector(new double[factDim]); double pred = predict(ins, orc, ufactSum, ifactSum); RealVector leftGrad = ifactSum; RealVector rightGrad = ufactSum; for (int i = 0; i < ins.ufeas.size(); i++) { orc.addVectorOracle(SVDFeatureKey.FACTORS.get(), ins.ufeas.get(i).getIndex(), leftGrad.mapMultiply(ins.ufeas.get(i).getValue())); }//w w w .ja v a2s .c o m for (int i = 0; i < ins.ifeas.size(); i++) { orc.addVectorOracle(SVDFeatureKey.FACTORS.get(), ins.ifeas.get(i).getIndex(), rightGrad.mapMultiply(ins.ifeas.get(i).getValue())); } orc.setValues(pred, ins.label, ins.weight); oracles.add(orc); } return oracles; }
From source file:org.lenskit.mf.BPR.BPRMFModelProvider.java
@Override public MFModel get() { // This will accumulate BPR-Opt (minus the regularization) and will be negated to make an error. // -30 is arbitrary, but would indicate a _really_ consistently bad prediction (~ p=1*10^-13), // and is therefore a reasonable "max_error". RollingWindowMeanAccumulator optAccum = new RollingWindowMeanAccumulator(10000, -30); // set up user index and matrix int userCount = dao.getEntityIds(CommonTypes.USER).size(); for (long uid : dao.getEntityIds(CommonTypes.USER)) { userIndex.internId(uid);/*ww w .j a va 2s. c om*/ } RealMatrix userFeatures = MatrixUtils.createRealMatrix(userCount, featureCount); initFeatures(userFeatures); // set up item index and matrix int itemCount = dao.getEntityIds(CommonTypes.ITEM).size(); for (long iid : dao.getEntityIds(CommonTypes.ITEM)) { itemIndex.internId(iid); } RealMatrix itemFeatures = MatrixUtils.createRealMatrix(itemCount, featureCount); initFeatures(itemFeatures); logger.debug("Learning rate is {}", learningRate); logger.debug("Regularization term is {}", regularization); logger.info("Building MPR-MF with {} features for {} users and {} items", featureCount, userCount, itemCount); TrainingLoopController controller = stoppingCondition.newLoop(); //REVIEW: because of the nature of training samples (and the point that the BPR paper makes that training // by-item or by-user are not optimal) one "iteration" here will be one training update. This leads to _really_ // big iteration counts, which can actually overflow ints!. one suggestion would be to allow the iteration count // controller to accept longs, but I don't know if this will introduce backwards compatibility issues (I imagine // this depends on the robustness of our type conversion in the configuration. while (controller.keepTraining(-optAccum.getMean())) { for (TrainingItemPair pair : pairGenerator.nextBatch()) { // Note: bad code style variable names are generally to match BPR paper and enable easier implementation long iid = pair.g; int i = itemIndex.internId(iid); long jid = pair.l; int j = itemIndex.internId(jid); long uid = pair.u; int u = userIndex.internId(uid); RealVector w_u = userFeatures.getRowVector(u); RealVector h_i = itemFeatures.getRowVector(i); RealVector h_j = itemFeatures.getRowVector(j); double xui = w_u.dotProduct(h_i); double xuj = w_u.dotProduct(h_j); double xuij = xui - xuj; double bprTerm = 1 / (1 + exp(xuij)); // w_u update RealVector h_i_j = h_i.subtract(h_j); RealVector w_u_update = w_u.mapMultiply(-regularization); w_u_update.combineToSelf(1, bprTerm, h_i_j); // h_i update RealVector h_i_update = h_i.mapMultiply(-regularization); h_i_update.combineToSelf(1, bprTerm, w_u); // h_j update RealVector h_j_update = h_j.mapMultiply(-regularization); h_j_update.combineToSelf(1, -bprTerm, w_u); // perform updates w_u.combineToSelf(1, learningRate, w_u_update); h_i.combineToSelf(1, learningRate, h_i_update); h_j.combineToSelf(1, learningRate, h_j_update); // update the optimization function accumulator (note we are not including the regularization term) optAccum.add(Math.log(1 / (1 + Math.exp(-xuij)))); } } return new MFModel(userFeatures, itemFeatures, userIndex, itemIndex); }
From source file:org.rhwlab.variationalbayesian.GaussianMixture.java
public void statistics() { // calculate N calculateN();// ww w.j a v a 2s . c om // calculate xBar for (int k = 0; k < K; ++k) { if (inModel[k]) { xBar[k].set(0.0); for (int n = 0; n < X.getN(); ++n) { Voxel vox = X.get(n); double intensity = vox.getAdjusted(); RealVector x = vox.coords; // double rnk = (intensity-background)*r.getEntry(n, k); double rnk = intensity * r.getEntry(n, k); xBar[k] = xBar[k].add(x.mapMultiply(rnk)); } xBar[k].mapDivideToSelf(N[k]); } } // calculate the S for (int k = 0; k < K; ++k) { if (inModel[k]) { if (k == 31) { int iusahdfis = 90; } for (int row = 0; row < X.getD(); ++row) { for (int col = 0; col < X.getD(); ++col) { S[k].setEntry(row, col, 0.0); } } for (int n = 0; n < X.getN(); ++n) { Voxel vox = X.get(n); RealVector x = vox.coords; RealVector del = x.subtract(xBar[k]); RealMatrix mat = del.outerProduct(del); double rnk = r.getEntry(n, k); if (rnk > 0) { int iuhsafuid = 0; } // S[k] = S[k].add(mat.scalarMultiply((intensity-background)*rnk)); S[k] = S[k].add(mat.scalarMultiply(vox.getAdjusted() * rnk)); } S[k] = S[k].scalarMultiply(1.0 / N[k]); } } int uisahdfius = 0; }
From source file:rcdemo.physics.FrictionForceModel.java
@Override public RealVector getForce(RealVector x, RealVector v) { double rho = 1; double Cd = 1; double A = 1; double factor = 0.5 * rho * v.getNorm() * Cd * A; return v.mapMultiply(-factor); }
From source file:rcdemo.track.TrackODE.java
@Override public void computeSecondDerivatives(double t, double[] y, double[] yDot, double[] yDDot) { double s = y[0]; double dsdt = yDot[0]; RealVector x = track.getx(s);//from ww w . j a v a 2 s . c om RealVector dxds = track.getDxDs(s); RealVector ddxdss = track.getDDxDss(s); RealVector u = dxds; RealVector v = dxds.mapMultiply(dsdt); RealVector F = forceModel.getForce(x, v); double ddsdtt = F.subtract(ddxdss.mapMultiply(dsdt * dsdt)).dotProduct(u) / dxds.dotProduct(u); //double E = 0.5 * v.dotProduct(v) + forceModel.getPotentialEnergy(x, v); //System.out.format("%4.1f: %s %s %s\n", t, s, x, E); //System.out.format("%7.5f: %s %s %s\n", t, s, x, E); yDDot[0] = ddsdtt; }