Examples with Array2DRowRealMatrix org.apache.commons.math3.linear.Array2DRowRealMatrix used on opensource projects

Example 61 with Array2DRowRealMatrix

use of org.apache.commons.math3.linear.Array2DRowRealMatrix in project narchy by automenta.

the class MyCMAESOptimizer method initializeCMA.

/**
 * Initialization of the dynamic search parameters
 *
 * @param guess Initial guess for the arguments of the fitness function.
 */
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);
    // overall standard deviation
    sigma = max(insigma);
    // initialize termination criteria
    stopTolUpX = oNEtHOUSAND * max(insigma);
    stopTolX = epsilonWTF11 * max(insigma);
    this.stopTolFun = EPSILON_WTF12;
    this.stopTolHistFun = epsilonwtf13;
    // initialize selection strategy parameters
    // number of parents/points for recombination
    mu = lambda / 2;
    /* log(mu + 0.5), stored for efficiency. */
    double logMu2 = Math.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);
    // variance-effectiveness of sum w_i x_i
    mueff = sumw * sumw / sumwq;
    // initialize dynamic strategy parameters and constants
    cc = (4 + mueff / dimension) / (dimension + 4 + 2 * mueff / dimension);
    cs = (mueff + 2) / (dimension + mueff + 3.0);
    damps = (1 + 2 * Math.max(0, Math.sqrt((mueff - 1) / (dimension + 1)) - 1)) * Math.max(0.3, 1 - dimension / (epsilon6WTF + maxIterations)) + // minor increment
    cs;
    ccov1 = 2 / ((dimension + 1.3) * (dimension + 1.3) + mueff);
    ccovmu = Math.min(1 - ccov1, 2 * (mueff - 2 + 1 / mueff) / ((dimension + 2) * (dimension + 2) + mueff));
    ccov1Sep = Math.min(1, ccov1 * (dimension + 1.5) / 3);
    ccovmuSep = Math.min(1 - ccov1, ccovmu * (dimension + 1.5) / 3);
    chiN = Math.sqrt(dimension) * (1 - 1 / (4.0 * dimension) + 1 / (21.0 * dimension * dimension));
    // intialize CMA internal values - updated each generation
    // objective variables
    xmean = MatrixUtils.createColumnRealMatrix(guess);
    diagD = insigma.scalarMultiply(1 / sigma);
    diagC = square(diagD);
    // evolution paths for C and sigma
    pc = zeros(dimension, 1);
    // B defines the coordinate system
    ps = zeros(dimension, 1);
    normps = ps.getFrobeniusNorm();
    B = eye(dimension, dimension);
    // diagonal D defines the scaling
    D = ones(dimension, 1);
    BD = times(B, repmat(diagD.transpose(), dimension, 1));
    // covariance
    C = B.multiply(diag(square(D)).multiply(B.transpose()));
    /* Size of history queue of best values. */
    int historySize = 10 + (int) (3 * 10 * dimension / (double) lambda);
    // history of fitness values
    fitnessHistory = new double[historySize];
    for (int i = 0; i < historySize; i++) {
        fitnessHistory[i] = Double.POSITIVE_INFINITY;
    }
}

Example 62 with Array2DRowRealMatrix

use of org.apache.commons.math3.linear.Array2DRowRealMatrix in project MindsEye by SimiaCryptus.

the class ObjectLocation method run.

/**
 * Run.
 *
 * @param log the log
 */
public void run(@Nonnull final NotebookOutput log) {
    @Nonnull String logName = "cuda_" + log.getName() + ".log";
    log.p(log.file((String) null, logName, "GPU Log"));
    CudaSystem.addLog(new PrintStream(log.file(logName)));
    ImageClassifier classifier = getClassifierNetwork();
    Layer classifyNetwork = classifier.getNetwork();
    ImageClassifier locator = getLocatorNetwork();
    Layer locatorNetwork = locator.getNetwork();
    ArtistryUtil.setPrecision((DAGNetwork) classifyNetwork, Precision.Float);
    ArtistryUtil.setPrecision((DAGNetwork) locatorNetwork, Precision.Float);
    Tensor[][] inputData = loadImages_library();
    // Tensor[][] inputData = loadImage_Caltech101(log);
    double alphaPower = 0.8;
    final AtomicInteger index = new AtomicInteger(0);
    Arrays.stream(inputData).limit(10).forEach(row -> {
        log.h3("Image " + index.getAndIncrement());
        final Tensor img = row[0];
        log.p(log.image(img.toImage(), ""));
        Result classifyResult = classifyNetwork.eval(new MutableResult(row));
        Result locationResult = locatorNetwork.eval(new MutableResult(row));
        Tensor classification = classifyResult.getData().get(0);
        List<CharSequence> categories = classifier.getCategories();
        int[] sortedIndices = IntStream.range(0, categories.size()).mapToObj(x -> x).sorted(Comparator.comparing(i -> -classification.get(i))).mapToInt(x -> x).limit(10).toArray();
        logger.info(Arrays.stream(sortedIndices).mapToObj(i -> String.format("%s: %s = %s%%", i, categories.get(i), classification.get(i) * 100)).reduce((a, b) -> a + "\n" + b).orElse(""));
        Map<CharSequence, Tensor> vectors = new HashMap<>();
        List<CharSequence> predictionList = Arrays.stream(sortedIndices).mapToObj(categories::get).collect(Collectors.toList());
        Arrays.stream(sortedIndices).limit(10).forEach(category -> {
            CharSequence name = categories.get(category);
            log.h3(name);
            Tensor alphaTensor = renderAlpha(alphaPower, img, locationResult, classification, category);
            log.p(log.image(img.toRgbImageAlphaMask(0, 1, 2, alphaTensor), ""));
            vectors.put(name, alphaTensor.unit());
        });
        Tensor avgDetection = vectors.values().stream().reduce((a, b) -> a.add(b)).get().scale(1.0 / vectors.size());
        Array2DRowRealMatrix covarianceMatrix = new Array2DRowRealMatrix(predictionList.size(), predictionList.size());
        for (int x = 0; x < predictionList.size(); x++) {
            for (int y = 0; y < predictionList.size(); y++) {
                Tensor l = vectors.get(predictionList.get(x)).minus(avgDetection);
                Tensor r = vectors.get(predictionList.get(y)).minus(avgDetection);
                covarianceMatrix.setEntry(x, y, l.dot(r));
            }
        }
        @Nonnull final EigenDecomposition decomposition = new EigenDecomposition(covarianceMatrix);
        for (int objectVector = 0; objectVector < 10; objectVector++) {
            log.h3("Eigenobject " + objectVector);
            double eigenvalue = decomposition.getRealEigenvalue(objectVector);
            RealVector eigenvector = decomposition.getEigenvector(objectVector);
            Tensor detectionRegion = IntStream.range(0, eigenvector.getDimension()).mapToObj(i -> vectors.get(predictionList.get(i)).scale(eigenvector.getEntry(i))).reduce((a, b) -> a.add(b)).get();
            detectionRegion = detectionRegion.scale(255.0 / detectionRegion.rms());
            CharSequence categorization = IntStream.range(0, eigenvector.getDimension()).mapToObj(i -> {
                CharSequence category = predictionList.get(i);
                double component = eigenvector.getEntry(i);
                return String.format("<li>%s = %.4f</li>", category, component);
            }).reduce((a, b) -> a + "" + b).get();
            log.p(String.format("Object Detected: <ol>%s</ol>", categorization));
            log.p("Object Eigenvalue: " + eigenvalue);
            log.p("Object Region: " + log.image(img.toRgbImageAlphaMask(0, 1, 2, detectionRegion), ""));
            log.p("Object Region Compliment: " + log.image(img.toRgbImageAlphaMask(0, 1, 2, detectionRegion.scale(-1)), ""));
        }
        // final int[] orderedVectors = IntStream.range(0, 10).mapToObj(x -> x)
        // .sorted(Comparator.comparing(x -> -decomposition.getRealEigenvalue(x))).mapToInt(x -> x).toArray();
        // IntStream.range(0, orderedVectors.length)
        // .mapToObj(i -> {
        // //double realEigenvalue = decomposition.getRealEigenvalue(orderedVectors[i]);
        // return decomposition.getEigenvector(orderedVectors[i]).toArray();
        // }
        // ).toArray(i -> new double[i][]);
        log.p(String.format("<table><tr><th>Cosine Distance</th>%s</tr>%s</table>", Arrays.stream(sortedIndices).limit(10).mapToObj(col -> "<th>" + categories.get(col) + "</th>").reduce((a, b) -> a + b).get(), Arrays.stream(sortedIndices).limit(10).mapToObj(r -> {
            return String.format("<tr><td>%s</td>%s</tr>", categories.get(r), Arrays.stream(sortedIndices).limit(10).mapToObj(col -> {
                return String.format("<td>%.4f</td>", Math.acos(vectors.get(categories.get(r)).dot(vectors.get(categories.get(col)))));
            }).reduce((a, b) -> a + b).get());
        }).reduce((a, b) -> a + b).orElse("")));
    });
    log.setFrontMatterProperty("status", "OK");
}

Example 63 with Array2DRowRealMatrix

use of org.apache.commons.math3.linear.Array2DRowRealMatrix in project MindsEye by SimiaCryptus.

the class ArtistryUtil method pca.

/**
 * Pca tensor.
 *
 * @param cov   the cov
 * @param power the power
 * @return the tensor
 */
@Nonnull
public static Tensor pca(final Tensor cov, final double power) {
    final int inputbands = (int) Math.sqrt(cov.getDimensions()[2]);
    final int outputbands = inputbands;
    Array2DRowRealMatrix realMatrix = new Array2DRowRealMatrix(inputbands, inputbands);
    cov.coordStream(false).forEach(c -> {
        double v = cov.get(c);
        int x = c.getIndex() % inputbands;
        int y = (c.getIndex() - x) / inputbands;
        realMatrix.setEntry(x, y, v);
    });
    Tensor[] features = PCAUtil.pcaFeatures(realMatrix, outputbands, new int[] { 1, 1, inputbands }, power);
    Tensor kernel = new Tensor(1, 1, inputbands * outputbands);
    PCAUtil.populatePCAKernel_1(kernel, features);
    return kernel;
}

Example 64 with Array2DRowRealMatrix

use of org.apache.commons.math3.linear.Array2DRowRealMatrix in project GDSC-SMLM by aherbert.

the class ApacheLvmFitter method computeFit.

@Override
public FitStatus computeFit(double[] y, final double[] fx, double[] a, double[] parametersVariance) {
    final int n = y.length;
    try {
        // Different convergence thresholds seem to have no effect on the resulting fit, only the
        // number of
        // iterations for convergence
        final double initialStepBoundFactor = 100;
        final double costRelativeTolerance = 1e-10;
        final double parRelativeTolerance = 1e-10;
        final double orthoTolerance = 1e-10;
        final double threshold = Precision.SAFE_MIN;
        // Extract the parameters to be fitted
        final double[] initialSolution = getInitialSolution(a);
        // TODO - Pass in more advanced stopping criteria.
        // Create the target and weight arrays
        final double[] yd = new double[n];
        // final double[] w = new double[n];
        for (int i = 0; i < n; i++) {
            yd[i] = y[i];
        // w[i] = 1;
        }
        final LevenbergMarquardtOptimizer optimizer = new LevenbergMarquardtOptimizer(initialStepBoundFactor, costRelativeTolerance, parRelativeTolerance, orthoTolerance, threshold);
        // @formatter:off
        final LeastSquaresBuilder builder = new LeastSquaresBuilder().maxEvaluations(Integer.MAX_VALUE).maxIterations(getMaxEvaluations()).start(initialSolution).target(yd);
        if (function instanceof ExtendedNonLinearFunction && ((ExtendedNonLinearFunction) function).canComputeValuesAndJacobian()) {
            // Compute together, or each individually
            builder.model(new ValueAndJacobianFunction() {

                final ExtendedNonLinearFunction fun = (ExtendedNonLinearFunction) function;

                @Override
                public Pair<RealVector, RealMatrix> value(RealVector point) {
                    final double[] p = point.toArray();
                    final org.apache.commons.lang3.tuple.Pair<double[], double[][]> result = fun.computeValuesAndJacobian(p);
                    return new Pair<>(new ArrayRealVector(result.getKey(), false), new Array2DRowRealMatrix(result.getValue(), false));
                }

                @Override
                public RealVector computeValue(double[] params) {
                    return new ArrayRealVector(fun.computeValues(params), false);
                }

                @Override
                public RealMatrix computeJacobian(double[] params) {
                    return new Array2DRowRealMatrix(fun.computeJacobian(params), false);
                }
            });
        } else {
            // Compute separately
            builder.model(new MultivariateVectorFunctionWrapper((NonLinearFunction) function, a, n), new MultivariateMatrixFunctionWrapper((NonLinearFunction) function, a, n));
        }
        final LeastSquaresProblem problem = builder.build();
        final Optimum optimum = optimizer.optimize(problem);
        final double[] parameters = optimum.getPoint().toArray();
        setSolution(a, parameters);
        iterations = optimum.getIterations();
        evaluations = optimum.getEvaluations();
        if (parametersVariance != null) {
            // Set up the Jacobian.
            final RealMatrix j = optimum.getJacobian();
            // Compute transpose(J)J.
            final RealMatrix jTj = j.transpose().multiply(j);
            final double[][] data = (jTj instanceof Array2DRowRealMatrix) ? ((Array2DRowRealMatrix) jTj).getDataRef() : jTj.getData();
            final FisherInformationMatrix m = new FisherInformationMatrix(data);
            setDeviations(parametersVariance, m);
        }
        // Compute function value
        if (fx != null) {
            final ValueFunction function = (ValueFunction) this.function;
            function.initialise0(a);
            function.forEach(new ValueProcedure() {

                int index;

                @Override
                public void execute(double value) {
                    fx[index++] = value;
                }
            });
        }
        // As this is unweighted then we can do this to get the sum of squared residuals
        // This is the same as optimum.getCost() * optimum.getCost(); The getCost() function
        // just computes the dot product anyway.
        value = optimum.getResiduals().dotProduct(optimum.getResiduals());
    } catch (final TooManyEvaluationsException ex) {
        return FitStatus.TOO_MANY_EVALUATIONS;
    } catch (final TooManyIterationsException ex) {
        return FitStatus.TOO_MANY_ITERATIONS;
    } catch (final ConvergenceException ex) {
        // Occurs when QR decomposition fails - mark as a singular non-linear model (no solution)
        return FitStatus.SINGULAR_NON_LINEAR_MODEL;
    } catch (final Exception ex) {
        // TODO - Find out the other exceptions from the fitter and add return values to match.
        return FitStatus.UNKNOWN;
    }
    return FitStatus.OK;
}

Example 65 with Array2DRowRealMatrix

use of org.apache.commons.math3.linear.Array2DRowRealMatrix in project GDSC-SMLM by aherbert.

the class MultivariateGaussianMixtureExpectationMaximizationTest method getColumnMeans.

/**
 * Gets the column means. This is done using the same method as the means in the Apache Commons
 * Math Covariance class.
 *
 * @param data the data
 * @return the column means
 */
private static double[] getColumnMeans(double[][] data) {
    final Array2DRowRealMatrix m = new Array2DRowRealMatrix(data);
    final Mean mean = new Mean();
    return IntStream.range(0, data[0].length).mapToDouble(i -> mean.evaluate(m.getColumn(i))).toArray();
}