List of usage examples for org.apache.commons.math3.linear RealVector add
public RealVector add(RealVector v) throws DimensionMismatchException
From source file:edu.stanford.cfuller.colocalization3d.correction.Correction.java
/** * Applies an existing correction to a single x-y position in the Image plane. * * @param x The x-position at which to apply the correction. * @param y The y-position at which to apply the correction. * @return A RealVector containing 3 elements-- the magnitude of the correction in the x, y, and z dimensions, in that order. *//*from w ww . j a va 2 s . c o m*/ public RealVector correctPosition(double x, double y) throws UnableToCorrectException { RealVector corrections = new ArrayRealVector(3, 0.0); RealVector distsToCentroids = this.getPositionsForCorrection().getColumnVector(0).mapSubtract(x) .mapToSelf(new Power(2)); distsToCentroids = distsToCentroids .add(this.getPositionsForCorrection().getColumnVector(1).mapSubtract(y).mapToSelf(new Power(2))); distsToCentroids.mapToSelf(new Sqrt()); RealVector distRatio = distsToCentroids.ebeDivide(this.getDistanceCutoffs()); RealVector distRatioBin = new ArrayRealVector(distRatio.getDimension(), 0.0); for (int i = 0; i < distRatio.getDimension(); i++) { if (distRatio.getEntry(i) <= 1) distRatioBin.setEntry(i, 1.0); } RealVector weights = distRatio.map(new Power(2.0)).mapMultiplyToSelf(-3).mapAddToSelf(1) .add(distRatio.map(new Power(3.0)).mapMultiplyToSelf(2)); weights = weights.ebeMultiply(distRatioBin); double sumWeights = 0; int countWeights = 0; for (int i = 0; i < weights.getDimension(); i++) { if (weights.getEntry(i) > 0) { sumWeights += weights.getEntry(i); countWeights++; } } if (countWeights == 0) { // this means there were no points in the correction dataset near the position being corrected. throw (new UnableToCorrectException( "Incomplete coverage in correction dataset at (x,y) = (" + x + ", " + y + ").")); } RealMatrix cX = new Array2DRowRealMatrix(countWeights, this.getCorrectionX().getColumnDimension()); RealMatrix cY = new Array2DRowRealMatrix(countWeights, this.getCorrectionX().getColumnDimension()); RealMatrix cZ = new Array2DRowRealMatrix(countWeights, this.getCorrectionX().getColumnDimension()); RealVector xVec = new ArrayRealVector(countWeights, 0.0); RealVector yVec = new ArrayRealVector(countWeights, 0.0); RealVector keptWeights = new ArrayRealVector(countWeights, 0.0); int keptCounter = 0; for (int i = 0; i < weights.getDimension(); i++) { if (weights.getEntry(i) > 0) { cX.setRowVector(keptCounter, this.getCorrectionX().getRowVector(i)); cY.setRowVector(keptCounter, this.getCorrectionY().getRowVector(i)); cZ.setRowVector(keptCounter, this.getCorrectionZ().getRowVector(i)); xVec.setEntry(keptCounter, x - this.getPositionsForCorrection().getEntry(i, 0)); yVec.setEntry(keptCounter, y - this.getPositionsForCorrection().getEntry(i, 1)); keptWeights.setEntry(keptCounter, weights.getEntry(i)); keptCounter++; } } double xCorr = 0; double yCorr = 0; double zCorr = 0; RealMatrix allCorrectionParameters = new Array2DRowRealMatrix(countWeights, numberOfCorrectionParameters); RealVector ones = new ArrayRealVector(countWeights, 1.0); allCorrectionParameters.setColumnVector(0, ones); allCorrectionParameters.setColumnVector(1, xVec); allCorrectionParameters.setColumnVector(2, yVec); allCorrectionParameters.setColumnVector(3, xVec.map(new Power(2))); allCorrectionParameters.setColumnVector(4, yVec.map(new Power(2))); allCorrectionParameters.setColumnVector(5, xVec.ebeMultiply(yVec)); for (int i = 0; i < countWeights; i++) { xCorr += allCorrectionParameters.getRowVector(i).dotProduct(cX.getRowVector(i)) * keptWeights.getEntry(i); yCorr += allCorrectionParameters.getRowVector(i).dotProduct(cY.getRowVector(i)) * keptWeights.getEntry(i); zCorr += allCorrectionParameters.getRowVector(i).dotProduct(cZ.getRowVector(i)) * keptWeights.getEntry(i); } xCorr /= sumWeights; yCorr /= sumWeights; zCorr /= sumWeights; corrections.setEntry(0, xCorr); corrections.setEntry(1, yCorr); corrections.setEntry(2, zCorr); return corrections; }
From source file:com.joptimizer.functions.SOCPLogarithmicBarrier.java
public double[] gradient(double[] X) { RealVector x = new ArrayRealVector(X); RealVector ret = new ArrayRealVector(dim); for (int i = 0; i < socpConstraintParametersList.size(); i++) { SOCPConstraintParameters param = socpConstraintParametersList.get(i); double t = this.buildT(param, x); RealVector u = this.buildU(param, x); double t2uu = t * t - u.dotProduct(u); RealMatrix Jacob = this.buildJ(param, x); int k = u.getDimension(); RealVector G = new ArrayRealVector(k + 1); G.setSubVector(0, u);/*from w w w .j a va2 s . c om*/ G.setEntry(k, -t); RealVector ret_i = Jacob.operate(G).mapMultiply((2. / t2uu)); ret = ret.add(ret_i); } return ret.toArray(); }
From source file:edu.stanford.cfuller.colocalization3d.correction.PositionCorrector.java
/** * Applies an existing correction to the positions of a set of objects, using a specified reference * and correction channel./*from w w w . j a v a2 s . c o m*/ * * @param imageObjects A Vector containing all the ImageObjects to be corrected. * @param c The Correction object to be used, which could have been generated with determineCorrection, or loaded from disk. * @return A RealVector with one entry per ImageObject, containing the corrected scalar distance between the object in the reference channel and the channel being corrected. */ public RealVector applyCorrection(Correction c, java.util.List<ImageObject> imageObjects) { int referenceChannel = this.parameters.getIntValueForKey(REF_CH_PARAM); int channelToCorrect = this.parameters.getIntValueForKey(CORR_CH_PARAM); RealVector diffs = new ArrayRealVector(imageObjects.size(), 0.0); RealVector averageVectorDiffs = null; for (int i = 0; i < imageObjects.size(); i++) { diffs.setEntry(i, imageObjects.get(i).getScalarDifferenceBetweenChannels(referenceChannel, channelToCorrect, this.pixelToDistanceConversions)); if (i == 0) { averageVectorDiffs = imageObjects.get(i) .getVectorDifferenceBetweenChannels(referenceChannel, channelToCorrect).copy(); } else { averageVectorDiffs = averageVectorDiffs.add(imageObjects.get(i) .getVectorDifferenceBetweenChannels(referenceChannel, channelToCorrect).map(new Abs())); } } averageVectorDiffs.mapDivideToSelf(imageObjects.size()); averageVectorDiffs = averageVectorDiffs.ebeMultiply(this.pixelToDistanceConversions); java.util.logging.Logger.getLogger("edu.stanford.cfuller.colocalization3d") .info("mean separation uncorrected = " + diffs.getL1Norm() / diffs.getDimension()); java.util.logging.Logger.getLogger("edu.stanford.cfuller.colocalization3d") .info("mean separation components uncorrected = " + averageVectorDiffs.toString()); RealVector newDiffs = new ArrayRealVector(imageObjects.size(), 0.0); if (this.parameters.getBooleanValueForKey("correct_images")) { for (int i = 0; i < imageObjects.size(); i++) { try { newDiffs.setEntry(i, this.correctSingleObjectVectorDifference(c, imageObjects.get(i), referenceChannel, channelToCorrect).getNorm()); imageObjects.get(i).setCorrectionSuccessful(true); } catch (UnableToCorrectException e) { newDiffs.setEntry(i, -1.0 * Double.MAX_VALUE); imageObjects.get(i).setCorrectionSuccessful(false); } } java.util.logging.Logger.getLogger("edu.stanford.cfuller.colocalization3d") .info("mean separation corrected = " + newDiffs.getL1Norm() / newDiffs.getDimension()); } else { newDiffs = diffs; } return newDiffs; }
From source file:edu.stanford.cfuller.imageanalysistools.fitting.NelderMeadMinimizer.java
/** * Runs the minimization of the specified function starting from an initial simplex. * @param f The ObjectiveFunction to be minimized. * @param initialSimplex The initial simplex to use to start optimization, as might be returned from {@link #generateInitialSimplex} * @return The parameters at the function minimum in the same order as specified for each point on the simplex. *///from w ww.j av a 2s .c o m public RealVector optimize(ObjectiveFunction f, RealMatrix initialSimplex) { RealMatrix currentSimplex = initialSimplex.copy(); double currTolVal = 1.0e6; RealVector values = new ArrayRealVector(initialSimplex.getRowDimension(), 0.0); RealVector centerOfMass = new ArrayRealVector(initialSimplex.getColumnDimension(), 0.0); boolean shouldEvaluate = false; long iterCounter = 0; while (Math.abs(currTolVal) > this.tol) { int maxIndex = 0; int minIndex = 0; double maxValue = -1.0 * Double.MAX_VALUE; double minValue = Double.MAX_VALUE; double secondMaxValue = -1.0 * Double.MAX_VALUE; centerOfMass.mapMultiplyToSelf(0.0); if (shouldEvaluate) { for (int i = 0; i < currentSimplex.getRowDimension(); i++) { RealVector currRow = currentSimplex.getRowVector(i); values.setEntry(i, f.evaluate(currRow)); } } for (int i = 0; i < currentSimplex.getRowDimension(); i++) { double currValue = values.getEntry(i); if (currValue < minValue) { minValue = currValue; minIndex = i; } if (currValue > maxValue) { secondMaxValue = maxValue; maxValue = currValue; maxIndex = i; } else if (currValue > secondMaxValue) { secondMaxValue = currValue; } } for (int i = 0; i < currentSimplex.getRowDimension(); i++) { if (i == maxIndex) continue; centerOfMass = centerOfMass.add(currentSimplex.getRowVector(i)); } centerOfMass.mapDivideToSelf(currentSimplex.getRowDimension() - 1); RealVector oldPoint = currentSimplex.getRowVector(maxIndex); RealVector newPoint = centerOfMass.subtract(oldPoint).mapMultiplyToSelf(a).add(centerOfMass); // newpoint = COM + a*(COM-oldpoint) double newValue = f.evaluate(newPoint); if (newValue < secondMaxValue) { // success if (newValue < minValue) { // best found so far //expansion RealVector expPoint = centerOfMass.subtract(oldPoint).mapMultiplyToSelf(g).add(centerOfMass); double expValue = f.evaluate(expPoint); if (expValue < newValue) { currentSimplex.setRowVector(maxIndex, expPoint); currTolVal = 2.0 * (expValue - maxValue) / (1.0e-20 + expValue + maxValue); values.setEntry(maxIndex, expValue); shouldEvaluate = false; continue; } } //reflection currentSimplex.setRowVector(maxIndex, newPoint); currTolVal = 2.0 * (newValue - maxValue) / (1.0e-20 + newValue + maxValue); values.setEntry(maxIndex, newValue); shouldEvaluate = false; continue; } //contraction RealVector conPoint = centerOfMass.subtract(oldPoint).mapMultiplyToSelf(r).add(centerOfMass); double conValue = f.evaluate(conPoint); if (conValue < maxValue) { currentSimplex.setRowVector(maxIndex, conPoint); currTolVal = 2.0 * (conValue - maxValue) / (1.0e-20 + conValue + maxValue); values.setEntry(maxIndex, conValue); shouldEvaluate = false; continue; } //reduction for (int i = 0; i < currentSimplex.getRowDimension(); i++) { if (i == minIndex) continue; currentSimplex.setRowVector(i, currentSimplex.getRowVector(i).subtract(currentSimplex.getRowVector(minIndex)) .mapMultiplyToSelf(s).add(currentSimplex.getRowVector(minIndex))); } double redValue = f.evaluate(currentSimplex.getRowVector(maxIndex)); currTolVal = 2.0 * (redValue - maxValue) / (1.0e-20 + redValue + maxValue); shouldEvaluate = true; if (iterCounter++ > 100000) { System.out.println("stalled? tol: " + currTolVal + " minValue: " + minValue); } } double minValue = Double.MAX_VALUE; RealVector minVector = null; for (int i = 0; i < currentSimplex.getRowDimension(); i++) { values.setEntry(i, f.evaluate(currentSimplex.getRowVector(i))); if (values.getEntry(i) < minValue) { minValue = values.getEntry(i); minVector = currentSimplex.getRowVector(i); } } return minVector; }
From source file:edu.utexas.cs.tactex.utilityestimation.UtilityEstimatorDefaultForConsumption.java
/** * Core method for estimating utility/*from w w w . j a v a 2 s. co m*/ */ 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; }
From source file:org.grouplens.lenskit.diffusion.Iterative.IterativeDiffusionItemScorer.java
/** * Score items by computing predicted ratings. * * @see ItemScoreAlgorithm#scoreItems(ItemItemModel, org.grouplens.lenskit.vectors.SparseVector, org.grouplens.lenskit.vectors.MutableSparseVector, NeighborhoodScorer) *///from w w w . j a v a 2 s .c o m @Override public void score(long user, @Nonnull MutableSparseVector scores) { UserHistory<? extends Event> history = dao.getEventsForUser(user, summarizer.eventTypeWanted()); if (history == null) { history = History.forUser(user); } SparseVector summary = summarizer.summarize(history); VectorTransformation transform = normalizer.makeTransformation(user, summary); MutableSparseVector normed = summary.mutableCopy(); transform.apply(normed); scores.clear(); int numItems = 1682; //algorithm.scoreItems(model, normed, scores, scorer); int num_updates = 300; double update_rate = 1; double threshold = 0.01; RealVector z_out = diffusionMatrix.preMultiply(VectorUtils.toRealVector(numItems, normed)); boolean updated = true; LongSortedSet known = normed.keySet(); int count_iter = 0; for (int i = 0; i < num_updates && updated; i++) { updated = false; RealVector temp = diffusionMatrix.preMultiply(z_out); temp.mapMultiplyToSelf(z_out.getNorm() / temp.getNorm()); RealVector temp_diff = z_out.add(temp.mapMultiplyToSelf(-1.0)); for (int j = 0; j < numItems; j++) { if (!known.contains((long) (j + 1))) { //if the rating is not one of the known ones if (Math.abs(temp_diff.getEntry(j)) > threshold) { // if difference is large enough, update updated = true; z_out.setEntry(j, (1.0 - update_rate) * z_out.getEntry(j) + update_rate * temp.getEntry(j)); } } } count_iter++; } System.out.println(count_iter); LongSortedSet testDomain = scores.keyDomain(); //fill up the score vector for (int i = 0; i < numItems; i++) { if (testDomain.contains((long) (i + 1))) { scores.set((long) (i + 1), z_out.getEntry(i)); } } // untransform the scores transform.unapply(scores); System.out.println(scores); }
From source file:org.lambda3.indra.composition.SumVectorComposer.java
@Override public RealVector compose(List<RealVector> vectors) { if (vectors.isEmpty()) { return null; } else if (vectors.size() == 1) { return vectors.get(0); } else {/*from w w w .j a v a 2 s . c o m*/ RealVector sum = vectors.get(0).add(vectors.get(1)); for (int i = 2; i < vectors.size(); i++) { sum = sum.add(vectors.get(i)); } return sum; } }
From source file:org.lambda3.indra.core.composition.SumVectorComposer.java
@Override public RealVector compose(List<RealVector> vectors) { logger.trace("Composing {} vectors", vectors.size()); if (vectors.isEmpty()) { return null; } else if (vectors.size() == 1) { return vectors.get(0); } else {/*from w w w .ja va 2 s . com*/ RealVector sum = vectors.get(0).add(vectors.get(1)); for (int i = 2; i < vectors.size(); i++) { sum = sum.add(vectors.get(i)); } return sum; } }
From source file:org.rhwlab.BHC.LogNode.java
public double alternative(List<RealVector> data) { int n = data.size(); double nP = n + beta; if (n > maxN) { return 0.0; }//from www.j av a2 s .c o m int D = data.get(0).getDimension(); RealVector xSum = new ArrayRealVector(D); // a vector of zeros RealMatrix XXT = new Array2DRowRealMatrix(D, D); for (RealVector v : data) { xSum = xSum.add(v); RealMatrix v2 = v.outerProduct(v); XXT = XXT.add(v2); } RealMatrix Sp = S.add(XXT); Sp = Sp.add(m.outerProduct(m).scalarMultiply(beta * n / (nP))); Sp = Sp.add(xSum.outerProduct(xSum).scalarMultiply(1.0 / (nP))); Sp = Sp.subtract(m.outerProduct(xSum).add(xSum.outerProduct(m)).scalarMultiply(beta / (nP))); LUDecomposition ed = new LUDecomposition(Sp); if (!ed.getSolver().isNonSingular()) { System.exit(10); } return Utils.eln(ed.getDeterminant()); }
From source file:org.rhwlab.BHC.LogNode.java
public double logMarginalLikelihood(List<RealVector> data) { int n = data.size(); if (n > maxN) { return 0.0; }/* w ww . j a va2s . co m*/ int D = data.get(0).getDimension(); double rP = beta + n; Double lnrP = Utils.eln(rP); double nuP = nu + n; RealMatrix C = new Array2DRowRealMatrix(D, D); RealVector X = new ArrayRealVector(D); // a vector of zeros for (RealVector v : data) { X = X.add(v); RealMatrix v2 = v.outerProduct(v); C = C.add(v2); } RealVector mP = (m.mapMultiply(beta).add(X)).mapDivide(rP); RealMatrix Sp = C.add(S); RealMatrix rmmP = mP.outerProduct(mP).scalarMultiply(rP); Sp = Sp.add(rmm).subtract(rmmP); LUDecomposition ed = new LUDecomposition(Sp); if (!ed.getSolver().isNonSingular()) { System.exit(10); } double logDetSp = Utils.eln(ed.getDeterminant()); // double logDetSp = alternative(data); Double ret = Utils.elnPow(logPi, -n * D / 2.0); ret = Utils.elnMult(ret, Utils.elnPow(lnBeta, D / 2.0)); ret = Utils.elnMult(ret, -Utils.elnPow(lnrP, D / 2.0)); ret = Utils.elnMult(ret, logdetSnu); ret = Utils.elnMult(ret, -Utils.elnPow(logDetSp, nuP / 2.0)); for (int i = 1; i <= D; ++i) { ret = Utils.elnMult(ret, Gamma.logGamma((nuP + 1 - i) / 2.0)); ret = Utils.elnMult(ret, -Gamma.logGamma((nu + 1 - i) / 2.0)); } if (ThreadedAlgorithm.time < 175) { ret = ret + data.size() * logRSDLikelihood(); } return ret; }