List of usage examples for org.apache.commons.math3.linear RealVector mapMultiplyToSelf
public RealVector mapMultiplyToSelf(double d)
From source file:edu.byu.nlp.stats.DirichletMLEOptimizableTest.java
/** * Computes a Newton-Raphson update in-place to alpha. *///from ww w. ja v a 2 s. co m private static RealVector newtonRaphsonUpdate(final double[][] data, double[] alpha) { // We'll compute the gold-standard value the "long" way (taking the inverse of the Hessian) RealMatrix hessian = new Array2DRowRealMatrix(alpha.length, alpha.length); for (int r = 0; r < hessian.getRowDimension(); r++) { for (int c = 0; c < hessian.getColumnDimension(); c++) { hessian.addToEntry(r, c, data.length * Gamma.trigamma(DoubleArrays.sum(alpha))); if (r == c) { hessian.addToEntry(r, c, -data.length * Gamma.trigamma(alpha[r])); } } } RealVector derivative = new ArrayRealVector(alpha.length); for (int k = 0; k < alpha.length; k++) { derivative.setEntry(k, data.length * (Gamma.digamma(DoubleArrays.sum(alpha)) - Gamma.digamma(alpha[k]))); for (double[] theta : data) { derivative.addToEntry(k, theta[k]); } } RealMatrix hessianInverse = new LUDecomposition(hessian).getSolver().getInverse(); RealVector negDiff = hessianInverse.preMultiply(derivative); negDiff.mapMultiplyToSelf(-1.0); RealVector expected = new ArrayRealVector(alpha, true); return expected.add(negDiff); }
From source file:edu.oregonstate.eecs.mcplan.ml.SequentialProjectionHashLearner.java
@Override public void run() { final RealMatrix cov_reg = MatrixUtils.createRealIdentityMatrix(X.getRowDimension()) .scalarMultiply(shrinkage);//from w ww . j a v a2 s . com for (int k = 0; k < K; ++k) { System.out.println("k = " + k); System.out.println("\tCovariance"); final RealMatrix cov = Xi_.multiply(Xi_.transpose()).add(cov_reg); // System.out.println( cov ); System.out.println("\tM"); final RealMatrix M = cov; // XL.multiply( Si_ ).multiply( XLt ).add( cov.scalarMultiply( eta ) ); // TODO: You only need the largest eigenvalue, so the full // decomposition is a ton of unnecessary work. System.out.println("\tM[" + M.getRowDimension() + "x" + M.getColumnDimension() + "]"); final EigenDecomposition evd = new EigenDecomposition(M); final RealVector w = evd.getEigenvector(0); w.mapMultiplyToSelf(1.0 / w.getNorm()); // if( Math.abs( w.getNorm() - 1.0 ) > 1e-2 ) { // System.out.println( "! Non-unit eigenvector: ||w|| = " + w.getNorm() ); // System.out.println( Math.abs( w.getNorm() - 1.0 ) ); // assert( false ); // } W.add(w); final RealMatrix w_wt = w.outerProduct(w); final RealMatrix S_tilde = XLt.multiply(w_wt).multiply(XL); T(S_tilde, Si_); Si_ = Si_.subtract(S_tilde.scalarMultiply(alpha)); Xi_ = Xi_.subtract(w_wt.multiply(Xi_)); } }
From source file:edu.stanford.cfuller.imageanalysistools.clustering.ObjectClustering.java
/** * Applies basic clustering to an Image with objects. * * This will use the long-range gaussian filtering approach to assign clusters; objects sufficiently near to each other will be smeared into a single object and assigned to the same cluster. * * @param input An Image mask labeled such that each object in the Image is assigned a unique nonzero greylevel value. These should start at 1 and be consecutive. * @param original The original image (not currently used... this is here to maintain the interface with a previous version that used this image) * @param gaussianFiltered The mask with a long range Gaussian filter applied (as from {@link #gaussianFilterMask}). This is an optional parameter; * input null to have this automatically generated by the method. This parameter is * chiefly useful to save computation time when running the clutering multiple times. * This will be modified, so if planning to reuse the Gaussian filtered image, pass in a copy. *///from ww w . j a v a2 s. c o m public static Image doBasicClustering(WritableImage input, Image original, Image gaussianFiltered) { RelabelFilter rlf = new RelabelFilter(); LabelFilter lf = new LabelFilter(); MaskFilter mf = new MaskFilter(); mf.setReferenceImage(input); Histogram h_individualCentromeres = new Histogram(input); WritableImage origCopy = null; if (gaussianFiltered == null) { origCopy = gaussianFilterMask(input).getWritableInstance(); } else { origCopy = gaussianFiltered.getWritableInstance(); } lf.apply(origCopy); WritableImage mapped = ImageFactory.createWritable(origCopy); Histogram h_mapped_0 = new Histogram(origCopy); //first, find the centroid of each cluster org.apache.commons.math3.linear.RealVector centroids_x = new ArrayRealVector(h_mapped_0.getMaxValue() + 1); org.apache.commons.math3.linear.RealVector centroids_y = new ArrayRealVector(h_mapped_0.getMaxValue() + 1); org.apache.commons.math3.linear.RealVector counts = new ArrayRealVector(h_mapped_0.getMaxValue() + 1); centroids_x.mapMultiplyToSelf(0.0); centroids_y.mapMultiplyToSelf(0.0); counts.mapMultiplyToSelf(0.0); for (ImageCoordinate i : origCopy) { if (origCopy.getValue(i) > 0) { int value = (int) origCopy.getValue(i); centroids_x.setEntry(value, centroids_x.getEntry(value) + i.get(ImageCoordinate.X)); centroids_y.setEntry(value, centroids_y.getEntry(value) + i.get(ImageCoordinate.Y)); counts.setEntry(value, counts.getEntry(value) + 1); } } for (int i = 0; i < counts.getDimension(); i++) { if (counts.getEntry(i) == 0) { counts.setEntry(i, 1); centroids_x.setEntry(i, -1 * origCopy.getDimensionSizes().get(ImageCoordinate.X)); centroids_y.setEntry(i, -1 * origCopy.getDimensionSizes().get(ImageCoordinate.Y)); } centroids_x.setEntry(i, centroids_x.getEntry(i) / counts.getEntry(i)); centroids_y.setEntry(i, centroids_y.getEntry(i) / counts.getEntry(i)); } for (ImageCoordinate i : origCopy) { if (mapped.getValue(i) > 0 || input.getValue(i) == 0) continue; double minDistance = Double.MAX_VALUE; int minIndex = 0; for (int j = 0; j < centroids_x.getDimension(); j++) { double dist = Math.hypot(centroids_x.getEntry(j) - i.get(ImageCoordinate.X), centroids_y.getEntry(j) - i.get(ImageCoordinate.Y)); if (dist < minDistance) { minDistance = dist; minIndex = j; } } mapped.setValue(i, minIndex); } int[] centromereAssignments = new int[h_individualCentromeres.getMaxValue() + 1]; java.util.Arrays.fill(centromereAssignments, 0); for (ImageCoordinate i : mapped) { if (input.getValue(i) > 0) { int value = (int) input.getValue(i); if (centromereAssignments[value] > 0) { mapped.setValue(i, centromereAssignments[value]); } else { centromereAssignments[value] = (int) mapped.getValue(i); } } } mf.apply(mapped); origCopy.copy(mapped); mf.setReferenceImage(origCopy); mf.apply(input); rlf.apply(input); rlf.apply(origCopy); return origCopy; }
From source file:edu.oregonstate.eecs.mcplan.ml.GaussianMixtureModel.java
private void init() { final int step = Math.max(1, n_ / k_); final double unif = 1.0 / k_; double acc = 0.0; final RandomPermutationIterator<double[]> r = new RandomPermutationIterator<double[]>(data_, rng_); final RandomPermutationIterator<double[]> rrepeat = new RandomPermutationIterator<double[]>(data_, r.permutation());/*from w w w .j av a 2 s .c o m*/ for (int i = 0; i < k_; ++i) { final RealVector mu = new ArrayRealVector(d_); for (int j = 0; j < step; ++j) { final double[] x = r.next(); final RealVector v = new ArrayRealVector(x); mu.combineToSelf(1.0, 1.0, v); } final double Zinv = 1.0 / step; mu.mapMultiplyToSelf(Zinv); RealMatrix Sigma = new Array2DRowRealMatrix(d_, d_); for (int j = 0; j < step; ++j) { final double[] x = rrepeat.next(); final RealVector v = new ArrayRealVector(x); v.combineToSelf(1.0, -1.0, mu); Sigma = Sigma.add(v.outerProduct(v)); } Sigma = Sigma.scalarMultiply(Zinv); pi_[i] = unif; acc += unif; mu_[i] = mu; Sigma_[i] = Sigma; //MatrixUtils.createRealIdentityMatrix( d_ ); } pi_[k_ - 1] += (1.0 - acc); // Round-off error }
From source file:edu.stanford.cfuller.imageanalysistools.filter.LocalBackgroundEstimationFilter.java
/** * Applies the LocalBackgroundEstimationFilter to an Image. * @param im The Image that will be replaced by the output Image. This can be anything of the correct dimensions except a shallow copy of the reference Image. *//*w w w.j a v a2 s . c o m*/ @Override public void apply(WritableImage im) { if (this.referenceImage == null) { throw new ReferenceImageRequiredException( "LocalBackgroundEstimationFilter requires a reference image."); } edu.stanford.cfuller.imageanalysistools.image.Histogram h = new edu.stanford.cfuller.imageanalysistools.image.Histogram( this.referenceImage); int numPixelsInBox = boxSize * boxSize; ImageCoordinate ic = this.referenceImage.getDimensionSizes(); ImageCoordinate icnew = ImageCoordinate.createCoordXYZCT(ic.get(ImageCoordinate.X) + 2 * boxSize, ic.get(ImageCoordinate.Y) + 2 * boxSize, ic.get(ImageCoordinate.Z), ic.get(ImageCoordinate.C), ic.get(ImageCoordinate.T)); WritableImage padded = ImageFactory.createWritable(icnew, -1.0f); float maxValue = 0; for (ImageCoordinate i : this.referenceImage) { icnew.quickSet(ImageCoordinate.X, i.quickGet(ImageCoordinate.X) + boxSize); icnew.quickSet(ImageCoordinate.Y, i.quickGet(ImageCoordinate.Y) + boxSize); icnew.quickSet(ImageCoordinate.Z, i.quickGet(ImageCoordinate.Z)); icnew.quickSet(ImageCoordinate.C, i.quickGet(ImageCoordinate.C)); icnew.quickSet(ImageCoordinate.T, i.quickGet(ImageCoordinate.T)); padded.setValue(icnew, this.referenceImage.getValue(i)); if (this.referenceImage.getValue(i) > maxValue) maxValue = this.referenceImage.getValue(i); } RealVector overallCounts = new org.apache.commons.math3.linear.ArrayRealVector(h.getMaxValue() + 1); RealMatrix countsByRow = new org.apache.commons.math3.linear.Array2DRowRealMatrix(2 * boxSize + 1, h.getMaxValue() + 1); //loop over columns for (int i = boxSize; i < im.getDimensionSizes().quickGet(ImageCoordinate.X) + boxSize; i++) { overallCounts.mapMultiplyToSelf(0.0); double[] overallCounts_a = overallCounts.toArray(); countsByRow = countsByRow.scalarMultiply(0.0); double[][] countsByRow_a = countsByRow.getData(); int countsByRow_rowZero_pointer = 0; for (int m = i - boxSize; m <= i + boxSize; m++) { for (int n = 0; n < 2 * boxSize + 1; n++) { icnew.quickSet(ImageCoordinate.X, m); icnew.quickSet(ImageCoordinate.Y, n); int value = (int) padded.getValue(icnew); if (value == -1) continue; overallCounts_a[value]++; countsByRow_a[(n + countsByRow_rowZero_pointer) % countsByRow.getRowDimension()][value]++; } } int currMedian = 0; int runningSum = 0; int k = 0; while (runningSum < (numPixelsInBox >> 1)) { runningSum += overallCounts_a[k++]; } runningSum -= overallCounts_a[k - 1]; currMedian = k - 1; icnew.quickSet(ImageCoordinate.X, i - boxSize); icnew.quickSet(ImageCoordinate.Y, 0); im.setValue(icnew, currMedian); int num_rows = countsByRow.getRowDimension(); for (int j = boxSize + 1; j < im.getDimensionSizes().quickGet(ImageCoordinate.Y) + boxSize; j++) { double[] toRemove = countsByRow_a[(countsByRow_rowZero_pointer) % num_rows]; for (int oc_counter = 0; oc_counter < overallCounts_a.length; oc_counter++) { overallCounts_a[oc_counter] -= toRemove[oc_counter]; } for (int c = 0; c < toRemove.length; c++) { if (c < currMedian) { runningSum -= toRemove[c]; } countsByRow_a[countsByRow_rowZero_pointer % num_rows][c] *= 0.0; } countsByRow_rowZero_pointer++; for (int c = i - boxSize; c <= i + boxSize; c++) { icnew.quickSet(ImageCoordinate.X, c); icnew.quickSet(ImageCoordinate.Y, j + boxSize); int value = (int) padded.getValue(icnew); if (value == -1) continue; countsByRow_a[(countsByRow_rowZero_pointer + num_rows - 1) % num_rows][value] += 1; overallCounts_a[value]++; if (value < currMedian) { runningSum++; } } //case 1: runningSum > half of box if (runningSum > (numPixelsInBox >> 1)) { k = currMedian - 1; while (runningSum > (numPixelsInBox >> 1)) { runningSum -= overallCounts_a[k--]; } currMedian = k + 1; } else if (runningSum < (numPixelsInBox >> 1)) { // case 2: runningSum < half of box k = currMedian; while (runningSum < (numPixelsInBox >> 1)) { runningSum += overallCounts_a[k++]; } currMedian = k - 1; runningSum -= overallCounts_a[k - 1]; } //cast 3: spot on, do nothing icnew.quickSet(ImageCoordinate.X, i - boxSize); icnew.quickSet(ImageCoordinate.Y, j - boxSize); im.setValue(icnew, currMedian); } } icnew.recycle(); }
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 ww w . j av a2 s . c om 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.stanford.cfuller.colocalization3d.correction.PositionCorrector.java
private RealVector correctSingleObjectVectorDifference(Correction c, ImageObject obj, int referenceChannel, int correctionChannel) throws UnableToCorrectException { RealVector corr = c.correctPosition(obj.getPositionForChannel(referenceChannel).getEntry(0), obj.getPositionForChannel(referenceChannel).getEntry(1)); boolean flip = this.parameters.getBooleanValueForKey("flip_channels_at_end"); if (flip)/*from w ww. java 2s .c o m*/ corr.mapMultiplyToSelf(-1.0); boolean invert_z = this.parameters.hasKeyAndTrue(INVERT_Z_PARAM); if (invert_z) corr.setEntry(2, -1.0 * corr.getEntry(2)); obj.applyCorrectionVectorToChannel(correctionChannel, corr); RealVector correctedVectorDiff = obj .getCorrectedVectorDifferenceBetweenChannels(referenceChannel, correctionChannel) .ebeMultiply(this.pixelToDistanceConversions); return correctedVectorDiff; }
From source file:edu.stanford.cfuller.imageanalysistools.fitting.DifferentialEvolutionMinimizer.java
/** * Performs a minimization of a function starting with a given population. * // w w w .ja v a2 s. c o m * @param population The population of parameters to start from, one population entry per row, one parameter per column. * @param f The function to be minimized. * @param parameterLowerBounds The lower bounds of each parameter. This must have the same size as the column dimension of the population. Generated parameter values less than these values will be discarded. * @param parameterUpperBounds The upper bounds of each paraemter. This must have the same size as the column dimension of the population. Generated parameter values greater than these values will be discarded. * @param populationSize The size of the population of parameters sets. This must be equal to the row dimension of the population. * @param scaleFactor Factor controlling the scale of crossed over entries during crossover events. * @param maxIterations The maximum number of iterations to allow before returning a result. * @param crossoverFrequency The frequency of crossover from 0 to 1. At any given parameter, this is the probability of initiating a crossover as well as the probability of ending one after it has started. * @param tol Relative function value tolerance controlling termination; if the maximal and minimal population values differ by less than this factor times the maximal value, optimization will terminate. * @return The parameter values at the minimal function value found. */ public RealVector minimizeWithPopulation(RealMatrix population, ObjectiveFunction f, RealVector parameterLowerBounds, RealVector parameterUpperBounds, int populationSize, double scaleFactor, int maxIterations, double crossoverFrequency, double tol) { int numberOfParameters = parameterLowerBounds.getDimension(); double currMinValue = Double.MAX_VALUE; double currMaxValue = -1.0 * Double.MAX_VALUE; int iterationCounter = maxIterations; double mutationProb = 0.01; // int totalIterations =0; RealVector values = new ArrayRealVector(populationSize); boolean[] valuesChanged = new boolean[populationSize]; java.util.Arrays.fill(valuesChanged, true); computeValues(f, population, values, valuesChanged); RealVector newVec = new ArrayRealVector(numberOfParameters); RealMatrix newPopulation = new Array2DRowRealMatrix(populationSize, numberOfParameters); while (iterationCounter > 0) { for (int i = 0; i < populationSize; i++) { int i1 = RandomGenerator.getGenerator().randInt(populationSize); int i2 = RandomGenerator.getGenerator().randInt(populationSize); int i3 = RandomGenerator.getGenerator().randInt(populationSize); newVec.mapMultiplyToSelf(0.0); boolean inBounds = true; boolean isCrossingOver = false; for (int j = 0; j < numberOfParameters; j++) { if ((RandomGenerator.rand() < crossoverFrequency) ^ isCrossingOver) { if (!isCrossingOver) { isCrossingOver = true; } newVec.setEntry(j, scaleFactor * (population.getEntry(i2, j) - population.getEntry(i1, j)) + population.getEntry(i3, j)); } else { if (isCrossingOver) { isCrossingOver = false; } newVec.setEntry(j, population.getEntry(i, j)); } //random 10% range +/- mutation if ((RandomGenerator.rand() < mutationProb)) { double magnitude = 0.2 * (parameterUpperBounds.getEntry(j) - parameterLowerBounds.getEntry(j)); newVec.setEntry(j, newVec.getEntry(j) + (RandomGenerator.rand() - 0.5) * magnitude); } if (newVec.getEntry(j) < parameterLowerBounds.getEntry(j) || newVec.getEntry(j) > parameterUpperBounds.getEntry(j)) { inBounds = false; } } double functionValue = Double.MAX_VALUE; if (inBounds) functionValue = f.evaluate(newVec); //if (inBounds) System.out.printf("in bounds candidate value: %1.2f old value: %1.2f \n", functionValue, values.getEntry(i)); if (functionValue < values.getEntry(i)) { newPopulation.setRowVector(i, newVec); valuesChanged[i] = true; values.setEntry(i, functionValue); } else { newPopulation.setRowVector(i, population.getRowVector(i)); valuesChanged[i] = false; } } population = newPopulation; double tempMinValue = Double.MAX_VALUE; double tempMaxValue = -1.0 * Double.MAX_VALUE; // double averageValue = 0; for (int i = 0; i < values.getDimension(); i++) { double value = values.getEntry(i); // averageValue += value; if (value < tempMinValue) { tempMinValue = value; } if (value > tempMaxValue) { tempMaxValue = value; } } // averageValue /= values.getDimension(); currMinValue = tempMinValue; currMaxValue = tempMaxValue; //LoggingUtilities.getLogger().info("Iteration counter: " + Integer.toString(totalIterations) + " best score: " + currMinValue + " worst score: " + currMaxValue + " average score: " + averageValue); if (Math.abs(currMaxValue - currMinValue) < Math.abs(tol * currMaxValue) + Math.abs(tol * currMinValue)) { iterationCounter--; } else { iterationCounter = 1; } // totalIterations++; } for (int i = 0; i < populationSize; i++) { valuesChanged[i] = true; } computeValues(f, population, values, valuesChanged); double tempMinValue = Double.MAX_VALUE; int tempMinIndex = 0; for (int i = 0; i < populationSize; i++) { if (values.getEntry(i) < tempMinValue) { tempMinValue = values.getEntry(i); tempMinIndex = i; } } RealVector output = new ArrayRealVector(numberOfParameters); for (int i = 0; i < numberOfParameters; i++) { output.setEntry(i, population.getEntry(tempMinIndex, i)); } return output; }
From source file:edu.oregonstate.eecs.mcplan.ml.GaussianMixtureModel.java
@Override public void run() { init();/*w w w . j av a2 s .c o m*/ System.out.println("Init"); for (int i = 0; i < mu().length; ++i) { System.out.println("Mu " + i + ": " + mu()[i]); System.out.println("Sigma " + i + ": " + Sigma()[i]); } int iterations = 0; while (!converged_ && iterations++ < max_iterations_) { // Expectation makeMixture(); for (int i = 0; i < n_; ++i) { for (int j = 0; j < k_; ++j) { c_[i][j] = posterior(data_[i], j); } Fn.normalize_inplace(c_[i]); } // Maximization for (int j = 0; j < k_; ++j) { double Z = 0.0; final RealVector mu_j = new ArrayRealVector(d_); RealMatrix Sigma_j = new Array2DRowRealMatrix(d_, d_); for (int i = 0; i < n_; ++i) { final double c_ij = c_[i][j]; Z += c_ij; final RealVector x_i = new ArrayRealVector(data_[i]); // mu_j += c_ij * x_i mu_j.combineToSelf(1.0, 1.0, x_i.mapMultiply(c_ij)); final RealVector v = x_i.subtract(mu_[j]); // Sigma_j += c_ij * |v><v| Sigma_j = Sigma_j.add(v.outerProduct(v).scalarMultiply(c_ij)); } final double Zinv = 1.0 / Z; final double pi_j = Z / n_; mu_j.mapMultiplyToSelf(Zinv); Sigma_j = Sigma_j.scalarMultiply(Zinv); // converged &= hasConverged( j, pi_j, mu_j, Sigma_j ); pi_[j] = pi_j; mu_[j] = mu_j; Sigma_[j] = Sigma_j; } // debug(); final double log_likelihood = logLikelihood(); if (Math.abs(log_likelihood - old_log_likelihood_) < epsilon_) { converged_ = true; } old_log_likelihood_ = log_likelihood; } }
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 ww w. j av a 2s . c om*/ @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); }