org.apache.commons.math3.linear.RealVector Java Examples

/**
 * Predict the internal state estimation one time step ahead.
 *
 * @param u
 *            the control vector
 * @throws DimensionMismatchException
 *             if the dimension of the control vector does not match
 */
public void predict(final RealVector u) throws DimensionMismatchException {
    // sanity checks
    if (u != null &&
        u.getDimension() != controlMatrix.getColumnDimension()) {
        throw new DimensionMismatchException(u.getDimension(),
                                             controlMatrix.getColumnDimension());
    }

    // project the state estimation ahead (a priori state)
    // xHat(k)- = A * xHat(k-1) + B * u(k-1)
    stateEstimation = transitionMatrix.operate(stateEstimation);

    // add control input if it is available
    if (u != null) {
        stateEstimation = stateEstimation.add(controlMatrix.operate(u));
    }

    // project the error covariance ahead
    // P(k)- = A * P(k-1) * A' + Q
    errorCovariance = transitionMatrix.multiply(errorCovariance)
            .multiply(transitionMatrixT)
            .add(processModel.getProcessNoise());
}
 

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()));
        }
        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;
}
 

/**
 * @return a comparator for sorting the optima.
 */
private Comparator<PointVectorValuePair> getPairComparator() {
    return new Comparator<PointVectorValuePair>() {
        private final RealVector target = new ArrayRealVector(optimizer.getTarget(), false);
        private final RealMatrix weight = optimizer.getWeight();

        public int compare(final PointVectorValuePair o1,
                           final PointVectorValuePair o2) {
            if (o1 == null) {
                return (o2 == null) ? 0 : 1;
            } else if (o2 == null) {
                return -1;
            }
            return Double.compare(weightedResidual(o1),
                                  weightedResidual(o2));
        }

        private double weightedResidual(final PointVectorValuePair pv) {
            final RealVector v = new ArrayRealVector(pv.getValueRef(), false);
            final RealVector r = target.subtract(v);
            return r.dotProduct(weight.operate(r));
        }
    };
}
 

private static void plot(RealVector vector, int height, int width, String imageFile) throws Exception{

        BufferedImage image = new BufferedImage(width,height,BufferedImage.TYPE_INT_RGB);
//        Graphics2D g2d = image.createGraphics();
//        g2d.setBackground(Color.WHITE);
//
//
//        g2d.fillRect ( 0, 0, image.getWidth(), image.getHeight() );
//        g2d.dispose();
        for (int i=0;i<width;i++){
            for (int j=0;j<height;j++){
                int v = (int)(vector.getEntry(i*width+j));
                int rgb = 65536 * v + 256 * v + v;
                image.setRGB(j,i,rgb);
//                image.setRGB(j,i,(int)(vector.get(i*width+j)/255*16777215));
            }
        }


        new File(imageFile).getParentFile().mkdirs();
        ImageIO.write(image,"png",new File(imageFile));
    }
 

@Test
public void nonZeroScoreTest() {
    List<MultivariateDistribution> listDist = new ArrayList<>(3);
    double[] weights = {2. / 7, 3. / 7, 2. / 7};
    double[][] distData = {
            {1.5, 2}, {0.5, 0.4, 0.4, 0.5}, {2000},
            {2, 0}, {0.3, 0, 0, 0.6}, {3000},
            {4.5, 1}, {0.9, 0.2, 0.2, 0.3}, {2000}};
    for (int i = 0; i < distData.length; i += 3) {
        RealVector mean = new ArrayRealVector(distData[i + 0]);
        double[][] covArray = new double[2][2];
        covArray[0] = Arrays.copyOfRange(distData[i + 1], 0, 2);
        covArray[1] = Arrays.copyOfRange(distData[i + 1], 2, 4);
        RealMatrix cov = new BlockRealMatrix(covArray);
        listDist.add(new MultivariateNormal(mean, cov));
    }

    Mixture mixture = new Mixture(listDist, weights);

    assertEquals(0.155359, mixture.density(new ArrayRealVector(distData[0])), 1e-6);
    assertEquals(0.162771, mixture.density(new ArrayRealVector(distData[3])), 1e-6);
    assertEquals(0.094819, mixture.density(new ArrayRealVector(distData[6])), 1e-6);
}
 

private Map<String, RealVector> collectVectors(Iterable<? extends String> terms, String targetField) {
    Map<String, RealVector> results = new LinkedHashMap<>();

    try {
        TopDocs topDocs = LuceneUtils.getTopDocs(termsSearcher, terms, targetField);

        if (topDocs != null && topDocs.totalHits > 0) {
            for (ScoreDoc sd : topDocs.scoreDocs) {
                Document doc = termsSearcher.doc(sd.doc);
                RealVector v = BinaryCodecs.unmarshall(doc.getBinaryValue(VECTOR_FIELD).bytes, true, (int) getMetadata().dimensions);
                results.put(doc.get(TERM_FIELD), v);
            }
        }
    } catch (IOException e) {
        logger.error(e.getMessage());
        //TODO throw new expection here.
    }

    return results;
}
 

Datum processRecord(Record rec) throws Exception {
    int idx = 0;
    RealVector metricVec = new ArrayRealVector(metrics.size());
    for (String metric : metrics) {
        metricVec.setEntry(idx, rec.values.get(metric));
        ++idx;
    }

    idx = 1;
    List<Integer> attrList = new ArrayList<>(attributes.size());
    for (String attr : attributes) {
        attrList.add(conf.getEncoder().getIntegerEncoding(idx, rec.attributes.get(attr)));
        ++idx;
    }

    return new Datum(attrList, metricVec);
}
 

/**
 * @return a comparator for sorting the optima.
 */
private Comparator<PointVectorValuePair> getPairComparator() {
    return new Comparator<PointVectorValuePair>() {
        private final RealVector target = new ArrayRealVector(optimizer.getTarget(), false);
        private final RealMatrix weight = optimizer.getWeight();

        public int compare(final PointVectorValuePair o1,
                           final PointVectorValuePair o2) {
            if (o1 == null) {
                return (o2 == null) ? 0 : 1;
            } else if (o2 == null) {
                return -1;
            }
            return Double.compare(weightedResidual(o1),
                                  weightedResidual(o2));
        }

        private double weightedResidual(final PointVectorValuePair pv) {
            final RealVector v = new ArrayRealVector(pv.getValueRef(), false);
            final RealVector r = target.subtract(v);
            return r.dotProduct(weight.operate(r));
        }
    };
}
 

/**
 * @param mu1 mean vector of the first normal distribution
 * @param sigma1 covariance matrix of the first normal distribution
 * @param mu2 mean vector of the second normal distribution
 * @param sigma2 covariance matrix of the second normal distribution
 * @return the Hellinger distance between two multivariate normal distributions
 * @link https://en.wikipedia.org/wiki/Hellinger_distance#Examples
 */
public static double hellingerDistance(@Nonnull final RealVector mu1,
        @Nonnull final RealMatrix sigma1, @Nonnull final RealVector mu2,
        @Nonnull final RealMatrix sigma2) {
    RealVector muSub = mu1.subtract(mu2);
    RealMatrix sigmaMean = sigma1.add(sigma2).scalarMultiply(0.5d);
    LUDecomposition LUsigmaMean = new LUDecomposition(sigmaMean);
    double denominator = Math.sqrt(LUsigmaMean.getDeterminant());
    if (denominator == 0.d) {
        return 1.d; // avoid divide by zero
    }
    RealMatrix sigmaMeanInv = LUsigmaMean.getSolver().getInverse(); // has inverse iff det != 0
    double sigma1Det = MatrixUtils.det(sigma1);
    double sigma2Det = MatrixUtils.det(sigma2);

    double numerator = Math.pow(sigma1Det, 0.25d) * Math.pow(sigma2Det, 0.25d)
            * Math.exp(-0.125d * sigmaMeanInv.preMultiply(muSub).dotProduct(muSub));
    return 1.d - numerator / denominator;
}
 

public double logDensity(RealVector instance){
    double[] arr = new double[numComponents];
    for (int k=0;k<numComponents;k++){
        arr[k]=Math.log(mixtureCoefficients[k])+gaussianDistributions[k].logDensity(instance);
    }
    return MathUtil.logSumExp(arr);
}
 

/**
 * Store the illumination and index of the brightest sample.
 * @param illuminationFromSample illumination received from sample
 * @param sample current sample illuminating the element
 */
void store(final double illuminationFromSample,
           final Map.Entry<RealVector, Double> sample) {
    if (illuminationFromSample > this.brightestIllumination) {
        this.brightestIllumination = illuminationFromSample;
        this.brightestSample = sample;
    }
}
 

public SumPaths(WeightedIntDiGraph wg, RealVector s, RealVector t, double lambda, double sigma) {
        this.lambda = lambda;
        this.sigma = sigma;
        startWeight = s;
        endWeight = t;
        this.wg = wg;
        Pair<IntDiGraph, IntObjectBimap<DiEdge>> p = wg.edgeGraph(false);
        edgeGraph = p.get1();
//        haltAction = edgeGraph.max() + 1;
        // TODO: hook the halt action in
//        edgeGraph.addNode(haltAction);
        this.edgesToNodes = p.get2();
        this.goldWeight = WeightedIntDiGraph.sumWalks(wg.toMatrix(), s, t);
    }
 

/**
 * Get the -1 times the sum of all coefficients in the given array.
 * @param coefficients coefficients to sum
 * @return the -1 times the sum of all coefficients in the given array.
 */
protected static double getInvertedCoefficientSum(final RealVector coefficients) {
    double sum = 0;
    for (double coefficient : coefficients.toArray()) {
        sum -= coefficient;
    }
    return sum;
}
 

/**
 * Verifies that newSampleData methods consistently insert unitary columns
 * in design matrix.  Confirms the fix for MATH-411.
 */
@Test
public void testNewSample() {
    double[] design = new double[] {
      1, 19, 22, 33,
      2, 20, 30, 40,
      3, 25, 35, 45,
      4, 27, 37, 47
    };
    double[] y = new double[] {1, 2, 3, 4}; 
    double[][] x = new double[][] {
      {19, 22, 33},
      {20, 30, 40},
      {25, 35, 45},
      {27, 37, 47}   
    };
    AbstractMultipleLinearRegression regression = createRegression();
    regression.newSampleData(design, 4, 3);
    RealMatrix flatX = regression.getX().copy();
    RealVector flatY = regression.getY().copy();
    regression.newXSampleData(x);
    regression.newYSampleData(y);
    Assert.assertEquals(flatX, regression.getX());
    Assert.assertEquals(flatY, regression.getY());
    
    // No intercept
    regression.setNoIntercept(true);
    regression.newSampleData(design, 4, 3);
    flatX = regression.getX().copy();
    flatY = regression.getY().copy();
    regression.newXSampleData(x);
    regression.newYSampleData(y);
    Assert.assertEquals(flatX, regression.getX());
    Assert.assertEquals(flatY, regression.getY());
}
 

@Override
public double density(RealVector u) {
    double rtn = 1.0;
    final int d = u.getDimension();
    for (int i = 0; i < d; i++) {
        double i2 = u.getEntry(i) * u.getEntry(i);
        if (i2 > 1) {
            return 0;
        }
        rtn *= 1 - i2;
    }
    return Math.pow(0.75, d) * rtn;
}
 

public void mStep(){
    IntStream.range(0,gmm.getNumComponents()).parallel()
            .forEach(k->{
                double sumGamma = computeSumGamma(k);
                gmm.setMixtureCoefficient(k,sumGamma/data.getRowDimension());
                RealVector mean = computeMean(k, sumGamma);
                gmm.getGaussianDistributions()[k].setMean(mean);
                RealMatrix cov = computeCov(k, mean, sumGamma);
                RealMatrix stabilizedCov = stablize(cov);
                gmm.getGaussianDistributions()[k].setCovariance(stabilizedCov);
            });
}
 

@Test
public void testGetIterations() {
    LeastSquaresProblem lsp = base()
            .target(new double[]{1})
            .weight(new DiagonalMatrix(new double[]{1}))
            .start(new double[]{3})
            .model(new MultivariateJacobianFunction() {
                public Pair<RealVector, RealMatrix> value(final RealVector point) {
                    return new Pair<RealVector, RealMatrix>(
                            new ArrayRealVector(
                                    new double[]{
                                            FastMath.pow(point.getEntry(0), 4)
                                    },
                                    false),
                            new Array2DRowRealMatrix(
                                    new double[][]{
                                            {0.25 * FastMath.pow(point.getEntry(0), 3)}
                                    },
                                    false)
                    );
                }
            })
            .build();

    Optimum optimum = optimizer.optimize(lsp);

    //TODO more specific test? could pass with 'return 1;'
    Assert.assertTrue(optimum.getIterations() > 0);
}
 

/**
 * Get the -1 times the sum of all coefficients in the given array.
 * @param coefficients coefficients to sum
 * @return the -1 times the sum of all coefficients in the given array.
 */
protected static double getInvertedCoefficientSum(final RealVector coefficients) {
    double sum = 0;
    for (double coefficient : coefficients.toArray()) {
        sum -= coefficient;
    }
    return sum;
}
 

/**
 * Calculates beta by GLS.
 * <pre>
 *  b=(X' Omega^-1 X)^-1X'Omega^-1 y
 * </pre>
 * @return beta
 */
@Override
protected RealVector calculateBeta() {
    RealMatrix OI = getOmegaInverse();
    RealMatrix XT = getX().transpose();
    RealMatrix XTOIX = XT.multiply(OI).multiply(getX());
    RealMatrix inverse = new LUDecomposition(XTOIX).getSolver().getInverse();
    return inverse.multiply(XT).multiply(OI).operate(getY());
}
 

@Override
protected RealVector solve(final RealMatrix jacobian,
                           final RealVector residuals) {
    try {
        final Pair<RealMatrix, RealVector> normalEquation =
                computeNormalMatrix(jacobian, residuals);
        final RealMatrix normal = normalEquation.getFirst();
        final RealVector jTr = normalEquation.getSecond();
        return new CholeskyDecomposition(
                normal, SINGULARITY_THRESHOLD, SINGULARITY_THRESHOLD)
                .getSolver()
                .solve(jTr);
    } catch (NonPositiveDefiniteMatrixException e) {
        throw new ConvergenceException(LocalizedFormats.UNABLE_TO_SOLVE_SINGULAR_PROBLEM, e);
    }
}
 

@Test
public void reduceToVectorKalmanFilterTest() throws Exception {
    MacroBaseConf conf = new MacroBaseConf()
            .set(MacroBaseConf.RANDOM_SEED, 4)
            .set(MacroBaseConf.DATA_LOADER_TYPE, "CSV_LOADER")
            .set(MacroBaseConf.CSV_COMPRESSION, CSVIngester.Compression.GZIP)
            .set(MacroBaseConf.CSV_INPUT_FILE, "src/test/resources/data/2gaussians-500points.csv.gz")
            .set(MacroBaseConf.METRICS, "XX, YY")
            .set(MacroBaseConf.ATTRIBUTES, "");
    List<Datum> data = Drainer.drainIngest(conf);
    assertEquals(500, data.size());

    KalmanVectorFilter f = new KalmanVectorFilter(new ArrayRealVector(2), 1e-6, 1);

    RealMatrix shapeMatrix = new BlockRealMatrix(2, 1);
    KalmanFlattenedMatrixFilter mf = new KalmanFlattenedMatrixFilter(shapeMatrix, 1e-6, 1);

    List<Datum> oneCluster = data.subList(201, 500);
    List<RealVector> vectorFiltered = oneCluster.stream().map(d -> f.step(d.metrics(), 1)).collect(Collectors.toList());
    List<RealMatrix> matrixFiltered = oneCluster.stream()
            .map(d -> mf.step(AlgebraUtils.reshapeMatrixByColumns(d.metrics(), shapeMatrix), 1))
            .collect(Collectors.toList());

    for (int i = 0; i < 10; i++) {
        int ri = conf.getRandom().nextInt(300);
        assertEquals(vectorFiltered.get(ri), AlgebraUtils.flattenMatrixByColumns(matrixFiltered.get(ri)));
    }
}
 

@Override
protected Map<String, RealVector> collectVectors(Iterable<? extends String> terms) {
    Map<String, RealVector> results = new HashMap<>();
    for (String term : terms) {
        if (localCache.keySet().contains(term)) {
            results.put(term, localCache.get(term));
        }
    }

    return results;
}
 

/**
 *
 * @param y     the response vector
 * @param x     the design matrix
 */
private RealMatrix computeBeta(RealVector y, RealMatrix x) {
    if (solver == Solver.QR) {
        return computeBetaQR(y, x);
    } else {
        final int n = x.getRowDimension();
        final int p = x.getColumnDimension();
        final int offset = hasIntercept() ? 1 : 0;
        final RealMatrix xT = x.transpose();
        final RealMatrix xTxInv = new LUDecomposition(xT.multiply(x)).getSolver().getInverse();
        final RealVector betaVector = xTxInv.multiply(xT).operate(y);
        final RealVector residuals = y.subtract(x.operate(betaVector));
        this.rss = residuals.dotProduct(residuals);
        this.errorVariance = rss / (n - p);
        this.stdError = Math.sqrt(errorVariance);
        this.residuals = createResidualsFrame(residuals);
        final RealMatrix covMatrix = xTxInv.scalarMultiply(errorVariance);
        final RealMatrix result = new Array2DRowRealMatrix(p, 2);
        if (hasIntercept()) {
            result.setEntry(0, 0, betaVector.getEntry(0));      //Intercept coefficient
            result.setEntry(0, 1, covMatrix.getEntry(0, 0));    //Intercept variance
        }
        for (int i = 0; i < getRegressors().size(); i++) {
            final int index = i + offset;
            final double variance = covMatrix.getEntry(index, index);
            result.setEntry(index, 1, variance);
            result.setEntry(index, 0, betaVector.getEntry(index));
        }
        return result;
    }
}
 

@Test
public void equivalent() {
    final String CORPUS = "simple";
    final String DENSE = "dense";
    final String SPARSE = "sparse";

    RawSpaceModel rsm = getRawSpaceModel(DENSE, CORPUS);

    try {
        VectorIterator iter = rsm.getVectorIterator();


        while (iter.hasNext()) {
            TermVector v = iter.next();
            if (!v.term.equalsIgnoreCase(StandardPreProcessorIterator.NUMBER_PLACEHOLDER)) {
                Assert.assertNotNull(v.content);
                RealVector rvo = RealVectorUtil.loosePrecision(v.content);

                Map<String, RealVector> results = driver.getVectors(createVectorRequest(DENSE, CORPUS, v.term));
                RealVector rvDense = results.get(v.term);

                results = driver.getVectors(createVectorRequest(SPARSE, CORPUS, v.term));
                RealVector rvSparse = results.get(v.term);

                Assert.assertEquals(rvDense.toArray(), rvSparse.toArray());
            }
        }
    } catch (IOException e) {
        e.printStackTrace();
        Assert.fail();
    }
}
 

/**
 * Verifies that newSampleData methods consistently insert unitary columns
 * in design matrix.  Confirms the fix for MATH-411.
 */
@Test
public void testNewSample() {
    double[] design = new double[] {
      1, 19, 22, 33,
      2, 20, 30, 40,
      3, 25, 35, 45,
      4, 27, 37, 47
    };
    double[] y = new double[] {1, 2, 3, 4}; 
    double[][] x = new double[][] {
      {19, 22, 33},
      {20, 30, 40},
      {25, 35, 45},
      {27, 37, 47}   
    };
    AbstractMultipleLinearRegression regression = createRegression();
    regression.newSampleData(design, 4, 3);
    RealMatrix flatX = regression.getX().copy();
    RealVector flatY = regression.getY().copy();
    regression.newXSampleData(x);
    regression.newYSampleData(y);
    Assert.assertEquals(flatX, regression.getX());
    Assert.assertEquals(flatY, regression.getY());
    
    // No intercept
    regression.setNoIntercept(true);
    regression.newSampleData(design, 4, 3);
    flatX = regression.getX().copy();
    flatY = regression.getY().copy();
    regression.newXSampleData(x);
    regression.newYSampleData(y);
    Assert.assertEquals(flatX, regression.getX());
    Assert.assertEquals(flatY, regression.getY());
}
 

@Test
public void nonExistingTerms() {
    final String NON_EXISTENT_TERM = "yyyyyyyywywywywy";
    List<String> terms = Arrays.asList(NON_EXISTENT_TERM, "amor");
    VectorRequest request = new VectorRequest().language("PT").terms(terms).termComposition("SUM").model("dense").
            mt(true).corpus("simple").translationComposition("AVERAGE");
    Map<String, RealVector> results = driver.getVectors(request);
    Assert.assertEquals(results.size(), terms.size());
    Assert.assertNull(results.get(NON_EXISTENT_TERM));
}
 

/**
 * Check the value of a vector.
 * @param tolerance the absolute tolerance of comparisons
 * @param actual the vector to test
 * @param expected the expected values
 */
public void assertEquals(double tolerance, RealVector actual, double... expected){
    for (int i = 0; i < expected.length; i++) {
        Assert.assertEquals(expected[i], actual.getEntry(i), tolerance);
    }
    Assert.assertEquals(expected.length, actual.getDimension());
}
 

/**
 * Learn multiplicative correction factors as a function of GC content.  Basically, learn a
 * regression curve of coverage vs. GC in order to divide by that curve later.
 *
 * @param readCounts           integer read counts for a single sample
 * @param intervalGCContent    GC content (assumed to be from 0.0 to 1.0) of genomic intervals for {@code readCounts}
 */
private GCBiasCorrector(final RealVector readCounts,
                        final double[] intervalGCContent) {
    Utils.nonNull(readCounts);
    Utils.nonNull(intervalGCContent);
    Utils.validateArg(intervalGCContent.length > 0, "Number of intervals must be positive.");
    Utils.validateArg(intervalGCContent.length == readCounts.getDimension(),
            "Number of intervals in read-counts matrix and GC-content array do not match.");

    final List<List<Double>> readCountsByGC = new ArrayList<>(NUMBER_OF_GC_BINS);
    IntStream.range(0, NUMBER_OF_GC_BINS).forEach(n -> readCountsByGC.add(new ArrayList<>()));
    IntStream.range(0, intervalGCContent.length).forEach(n -> readCountsByGC.get(gcContentToBinIndex(intervalGCContent[n])).add(readCounts.getEntry(n)));
    gcCorrectionFactors = calculateCorrectionFactors(readCountsByGC);
}
 

private static int getVoxConstrDOFNum(RealVector coeff){
    //return what dof is constrained if this is a vox constr; else return -1
    int ans = -1;
    for(int dof=0; dof<coeff.getDimension(); dof++){
        double val = coeff.getEntry(dof);
        if(Math.abs(val)>1e-10){//nonzero
            if(ans!=-1)//this is not the first constrained dof.  Must not be a voxel constraint.  
                return -1;
            if(Math.abs(val-1)>1e-10)//coeff 1 expected for vox constr
                return -1;
            ans = dof;
        }
    }
    return ans;
}
 

/** {@inheritDoc} */
public Pair<RealVector, RealMatrix> value(final RealVector point) {
    //TODO get array from RealVector without copying?
    final double[] p = point.toArray();

    // Evaluate.
    return new Pair<RealVector, RealMatrix>(computeValue(p),
                                            computeJacobian(p));
}