Examples with DoubleMatrix1D cern.colt.matrix.DoubleMatrix1D used on opensource projects

Example 16 with DoubleMatrix1D

use of cern.colt.matrix.DoubleMatrix1D in project tdq-studio-se by Talend.

the class SmpBlas method dgemv.

public void dgemv(final boolean transposeA, final double alpha, DoubleMatrix2D A, final DoubleMatrix1D x, final double beta, DoubleMatrix1D y) {
    /*
	split A, as follows:
	
			x x
			x
			x
	A
	xxx     x y
	xxx     x
	---     -
	xxx     x
	xxx     x
	---     -
	xxx     x

	*/
    if (transposeA) {
        dgemv(false, alpha, A.viewDice(), x, beta, y);
        return;
    }
    int m = A.rows();
    int n = A.columns();
    long flops = 2L * m * n;
    // each thread should process at least 30000 flops
    int noOfTasks = (int) Math.min(flops / 30000, this.maxThreads);
    int width = A.rows();
    noOfTasks = Math.min(width, noOfTasks);
    if (noOfTasks < 2) {
        // parallelization doesn't pay off (too much start up overhead)
        seqBlas.dgemv(transposeA, alpha, A, x, beta, y);
        return;
    }
    // set up concurrent tasks
    int span = width / noOfTasks;
    final FJTask[] subTasks = new FJTask[noOfTasks];
    for (int i = 0; i < noOfTasks; i++) {
        final int offset = i * span;
        // last span may be a bit larger
        if (i == noOfTasks - 1)
            span = width - span * i;
        // split A along rows into blocks
        final DoubleMatrix2D AA = A.viewPart(offset, 0, span, n);
        final DoubleMatrix1D yy = y.viewPart(offset, span);
        subTasks[i] = new FJTask() {

            public void run() {
                seqBlas.dgemv(transposeA, alpha, AA, x, beta, yy);
            // System.out.println("Hello "+offset);
            }
        };
    }
    // run tasks and wait for completion
    try {
        this.smp.taskGroup.invoke(new FJTask() {

            public void run() {
                coInvoke(subTasks);
            }
        });
    } catch (InterruptedException exc) {
    }
}

Example 17 with DoubleMatrix1D

use of cern.colt.matrix.DoubleMatrix1D in project tdq-studio-se by Talend.

the class Sorting method sort.

/**
 *Sorts the matrix rows into ascending order, according to the <i>natural ordering</i> of the matrix values in the given column.
 *The returned view is backed by this matrix, so changes in the returned view are reflected in this matrix, and vice-versa.
 *To sort ranges use sub-ranging views. To sort columns by rows, use dice views. To sort descending, use flip views ...
 *<p>
 *<b>Example:</b>
 *<table border="1" cellspacing="0">
 *  <tr nowrap>
 *	<td valign="top"><tt>4 x 2 matrix: <br>
 *	  7, 6<br>
 *	  5, 4<br>
 *	  3, 2<br>
 *	  1, 0 <br>
 *	  </tt></td>
 *	<td align="left" valign="top">
 *	  <p><tt>column = 0;<br>
 *		view = quickSort(matrix,column);<br>
 *		System.out.println(view); </tt><tt><br>
 *		==> </tt></p>
 *	  </td>
 *	<td valign="top">
 *	  <p><tt>4 x 2 matrix:<br>
 *		1, 0<br>
 *		3, 2<br>
 *		5, 4<br>
 *		7, 6</tt><br>
 *		The matrix IS NOT SORTED.<br>
 *		The new VIEW IS SORTED.</p>
 *	  </td>
 *  </tr>
 *</table>
 *
 *@param matrix the matrix to be sorted.
 *@param column the index of the column inducing the order.
 *@return a new matrix view having rows sorted by the given column.
 *		<b>Note that the original matrix is left unaffected.</b>
 *@throws IndexOutOfBoundsException if <tt>column < 0 || column >= matrix.columns()</tt>.
 */
public DoubleMatrix2D sort(DoubleMatrix2D matrix, int column) {
    if (column < 0 || column >= matrix.columns())
        throw new IndexOutOfBoundsException("column=" + column + ", matrix=" + Formatter.shape(matrix));
    // row indexes to reorder instead of matrix itself
    int[] rowIndexes = new int[matrix.rows()];
    for (int i = rowIndexes.length; --i >= 0; ) rowIndexes[i] = i;
    final DoubleMatrix1D col = matrix.viewColumn(column);
    IntComparator comp = new IntComparator() {

        public int compare(int a, int b) {
            double av = col.getQuick(a);
            double bv = col.getQuick(b);
            // swap NaNs to the end
            if (av != av || bv != bv)
                return compareNaN(av, bv);
            return av < bv ? -1 : (av == bv ? 0 : 1);
        }
    };
    runSort(rowIndexes, 0, rowIndexes.length, comp);
    // take all columns in the original order
    return matrix.viewSelection(rowIndexes, null);
}

Example 18 with DoubleMatrix1D

use of cern.colt.matrix.DoubleMatrix1D in project tdq-studio-se by Talend.

the class Sorting method zdemo5.

/**
 * Demonstrates sorting with precomputation of aggregates (median and sum of logarithms).
 */
public static void zdemo5(int rows, int columns, boolean print) {
    Sorting sort = quickSort;
    // for reliable benchmarks, call this method twice: once with small dummy parameters to "warm up" the jitter, then with your real work-load
    System.out.println("\n\n");
    System.out.print("now initializing... ");
    cern.colt.Timer timer = new cern.colt.Timer().start();
    final cern.jet.math.Functions F = cern.jet.math.Functions.functions;
    DoubleMatrix2D A = cern.colt.matrix.DoubleFactory2D.dense.make(rows, columns);
    // initialize randomly
    A.assign(new cern.jet.random.engine.DRand());
    timer.stop().display();
    // also benchmark copying in its several implementation flavours
    DoubleMatrix2D B = A.like();
    timer.reset().start();
    System.out.print("now copying... ");
    B.assign(A);
    timer.stop().display();
    timer.reset().start();
    System.out.print("now copying subrange... ");
    B.viewPart(0, 0, rows, columns).assign(A.viewPart(0, 0, rows, columns));
    timer.stop().display();
    // System.out.println(A);
    timer.reset().start();
    System.out.print("now copying selected... ");
    B.viewSelection(null, null).assign(A.viewSelection(null, null));
    timer.stop().display();
    System.out.print("now sorting - quick version with precomputation... ");
    timer.reset().start();
    // THE QUICK VERSION (takes some 10 secs)
    A = sort.sort(A, hep.aida.bin.BinFunctions1D.median);
    // A = sort.sort(A,hep.aida.bin.BinFunctions1D.sumLog);
    timer.stop().display();
    // so we just show the first 5 rows
    if (print) {
        int r = Math.min(rows, 5);
        hep.aida.bin.BinFunction1D[] funs = { hep.aida.bin.BinFunctions1D.median, hep.aida.bin.BinFunctions1D.sumLog, hep.aida.bin.BinFunctions1D.geometricMean };
        String[] rowNames = new String[r];
        String[] columnNames = new String[columns];
        for (int i = columns; --i >= 0; ) columnNames[i] = Integer.toString(i);
        for (int i = r; --i >= 0; ) rowNames[i] = Integer.toString(i);
        System.out.println("first part of sorted result = \n" + new cern.colt.matrix.doublealgo.Formatter("%G").toTitleString(A.viewPart(0, 0, r, columns), rowNames, columnNames, null, null, null, funs));
    }
    System.out.print("now sorting - slow version... ");
    A = B;
    cern.colt.matrix.doublealgo.DoubleMatrix1DComparator fun = new cern.colt.matrix.doublealgo.DoubleMatrix1DComparator() {

        public int compare(DoubleMatrix1D x, DoubleMatrix1D y) {
            double a = cern.colt.matrix.doublealgo.Statistic.bin(x).median();
            double b = cern.colt.matrix.doublealgo.Statistic.bin(y).median();
            // double b = y.aggregate(F.plus,F.log);
            return a < b ? -1 : (a == b) ? 0 : 1;
        }
    };
    timer.reset().start();
    A = sort.sort(A, fun);
    timer.stop().display();
}

Example 19 with DoubleMatrix1D

use of cern.colt.matrix.DoubleMatrix1D in project tdq-studio-se by Talend.

the class Sorting method sort.

/**
 *Sorts the matrix rows according to the order induced by the specified comparator.
 *The returned view is backed by this matrix, so changes in the returned view are reflected in this matrix, and vice-versa.
 *The algorithm compares two rows (1-d matrices) at a time, determinining whether one is smaller, equal or larger than the other.
 *To sort ranges use sub-ranging views. To sort columns by rows, use dice views. To sort descending, use flip views ...
 *<p>
 *<b>Example:</b>
 *<pre>
 *// sort by sum of values in a row
 *DoubleMatrix1DComparator comp = new DoubleMatrix1DComparator() {
 *&nbsp;&nbsp;&nbsp;public int compare(DoubleMatrix1D a, DoubleMatrix1D b) {
 *&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;double as = a.zSum(); double bs = b.zSum();
 *&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;return as < bs ? -1 : as == bs ? 0 : 1;
 *&nbsp;&nbsp;&nbsp;}
 *};
 *sorted = quickSort(matrix,comp);
 *</pre>
 *
 *@param matrix the matrix to be sorted.
 *@param c the comparator to determine the order.
 *@return a new matrix view having rows sorted as specified.
 *		<b>Note that the original matrix is left unaffected.</b>
 */
public DoubleMatrix2D sort(final DoubleMatrix2D matrix, final DoubleMatrix1DComparator c) {
    // row indexes to reorder instead of matrix itself
    int[] rowIndexes = new int[matrix.rows()];
    for (int i = rowIndexes.length; --i >= 0; ) rowIndexes[i] = i;
    // precompute views for speed
    final DoubleMatrix1D[] views = new DoubleMatrix1D[matrix.rows()];
    for (int i = views.length; --i >= 0; ) views[i] = matrix.viewRow(i);
    IntComparator comp = new IntComparator() {

        public int compare(int a, int b) {
            // return c.compare(matrix.viewRow(a), matrix.viewRow(b));
            return c.compare(views[a], views[b]);
        }
    };
    runSort(rowIndexes, 0, rowIndexes.length, comp);
    // take all columns in the original order
    return matrix.viewSelection(rowIndexes, null);
}

Example 20 with DoubleMatrix1D

use of cern.colt.matrix.DoubleMatrix1D in project Gemma by PavlidisLab.

the class LinearModelAnalyzer method doAnalysis.

/**
 * @param bioAssaySet       source data, could be a SubSet
 * @param dmatrix           data
 * @param samplesUsed       analyzed
 * @param factors           included in the model
 * @param subsetFactorValue null unless analyzing a subset (only used for book-keeping)
 * @return analysis, or null if there was a problem.
 */
private DifferentialExpressionAnalysis doAnalysis(BioAssaySet bioAssaySet, DifferentialExpressionAnalysisConfig config, ExpressionDataDoubleMatrix dmatrix, List<BioMaterial> samplesUsed, List<ExperimentalFactor> factors, FactorValue subsetFactorValue) {
    // We may want to change this to fall back to running normally, though the real fix is to just finish the ebayes implementation.
    if (config.getModerateStatistics() && dmatrix.hasMissingValues()) {
        throw new UnsupportedOperationException("Ebayes cannot be used when there are values missing in the data");
    }
    if (factors.isEmpty()) {
        LinearModelAnalyzer.log.error("Must provide at least one factor");
        return null;
    }
    if (samplesUsed.size() <= factors.size()) {
        LinearModelAnalyzer.log.error("Must have more samples than factors");
        return null;
    }
    final Map<String, Collection<ExperimentalFactor>> label2Factors = this.getRNames(factors);
    Map<ExperimentalFactor, FactorValue> baselineConditions = ExperimentalDesignUtils.getBaselineConditions(samplesUsed, factors);
    this.dropIncompleteFactors(samplesUsed, factors, baselineConditions);
    if (factors.isEmpty()) {
        LinearModelAnalyzer.log.error("Must provide at least one factor; they were all removed due to incomplete values");
        return null;
    }
    QuantitationType quantitationType = dmatrix.getQuantitationTypes().iterator().next();
    ExperimentalFactor interceptFactor = this.determineInterceptFactor(factors, quantitationType);
    /*
         * Build our factor terms, with interactions handled specially
         */
    List<String[]> interactionFactorLists = new ArrayList<>();
    ObjectMatrix<String, String, Object> designMatrix = ExperimentalDesignUtils.buildDesignMatrix(factors, samplesUsed, baselineConditions);
    config.setBaseLineFactorValues(baselineConditions);
    boolean oneSampleTTest = interceptFactor != null && factors.size() == 1;
    if (!oneSampleTTest) {
        this.buildModelFormula(config, label2Factors, interactionFactorLists);
    }
    // calculate library size before log transformation (FIXME we compute this twice)
    DoubleMatrix1D librarySize = MatrixStats.colSums(dmatrix.getMatrix());
    /*
         * FIXME: remove columns that are marked as outliers.
         */
    dmatrix = ExpressionDataDoubleMatrixUtil.filterAndLog2Transform(quantitationType, dmatrix);
    DoubleMatrix<CompositeSequence, BioMaterial> namedMatrix = dmatrix.getMatrix();
    if (LinearModelAnalyzer.log.isDebugEnabled())
        this.outputForDebugging(dmatrix, designMatrix);
    /*
         * PREPARATION FOR 'NATIVE' FITTING
         */
    DoubleMatrix<String, String> sNamedMatrix = this.makeDataMatrix(designMatrix, namedMatrix);
    DesignMatrix properDesignMatrix = this.makeDesignMatrix(designMatrix, interactionFactorLists, baselineConditions);
    /*
         * Run the analysis NOTE this can be simplified if we strip out R code.
         */
    final Map<String, LinearModelSummary> rawResults = this.runAnalysis(namedMatrix, sNamedMatrix, properDesignMatrix, librarySize, config);
    if (rawResults.size() == 0) {
        LinearModelAnalyzer.log.error("Got no results from the analysis");
        return null;
    }
    /*
         * Initialize data structures we need to hold results.
         */
    // this used to be a Set, but a List is much faster.
    Map<String, List<DifferentialExpressionAnalysisResult>> resultLists = new HashMap<>();
    Map<String, List<Double>> pvaluesForQvalue = new HashMap<>();
    for (String factorName : label2Factors.keySet()) {
        resultLists.put(factorName, new ArrayList<DifferentialExpressionAnalysisResult>());
        pvaluesForQvalue.put(factorName, new ArrayList<Double>());
    }
    for (String[] fs : interactionFactorLists) {
        String intF = StringUtils.join(fs, ":");
        resultLists.put(intF, new ArrayList<DifferentialExpressionAnalysisResult>());
        pvaluesForQvalue.put(intF, new ArrayList<Double>());
    }
    if (pvaluesForQvalue.isEmpty()) {
        LinearModelAnalyzer.log.warn("No results were obtained for the current stage of analysis.");
        return null;
    }
    /*
         * Create result objects for each model fit. Keeping things in order is important.
         */
    final Transformer rowNameExtractor = TransformerUtils.invokerTransformer("getId");
    boolean warned = false;
    int notUsable = 0;
    int processed = 0;
    for (CompositeSequence el : namedMatrix.getRowNames()) {
        if (++processed % 15000 == 0) {
            LinearModelAnalyzer.log.info("Processed results for " + processed + " elements ...");
        }
        LinearModelSummary lm = rawResults.get(rowNameExtractor.transform(el).toString());
        if (LinearModelAnalyzer.log.isDebugEnabled())
            LinearModelAnalyzer.log.debug(el.getName() + "\n" + lm);
        if (lm == null) {
            if (!warned) {
                LinearModelAnalyzer.log.warn("No result for " + el + ", further warnings suppressed");
                warned = true;
            }
            notUsable++;
            continue;
        }
        for (String factorName : label2Factors.keySet()) {
            if (!pvaluesForQvalue.containsKey(factorName)) {
                // was dropped.
                continue;
            }
            Double overallPValue;
            DifferentialExpressionAnalysisResult probeAnalysisResult = DifferentialExpressionAnalysisResult.Factory.newInstance();
            probeAnalysisResult.setProbe(el);
            if (lm.getCoefficients() == null) {
                // probeAnalysisResult.setPvalue( null );
                // pvaluesForQvalue.get( factorName ).add( overallPValue );
                // resultLists.get( factorName ).add( probeAnalysisResult );
                notUsable++;
                continue;
            }
            Collection<ExperimentalFactor> factorsForName = label2Factors.get(factorName);
            if (factorsForName.size() > 1) {
                /*
                     * Interactions
                     */
                if (factorsForName.size() > 2) {
                    LinearModelAnalyzer.log.error("Handling more than two-way interactions is not implemented");
                    return null;
                }
                assert factorName.contains(":");
                String[] factorNames = StringUtils.split(factorName, ":");
                assert factorNames.length == factorsForName.size();
                overallPValue = lm.getInteractionEffectP(factorNames);
                if (overallPValue != null && !Double.isNaN(overallPValue)) {
                    Map<String, Double> interactionContrastTStats = lm.getContrastTStats(factorName);
                    Map<String, Double> interactionContrastCoeffs = lm.getContrastCoefficients(factorName);
                    Map<String, Double> interactionContrastPValues = lm.getContrastPValues(factorName);
                    for (String term : interactionContrastPValues.keySet()) {
                        Double contrastPvalue = interactionContrastPValues.get(term);
                        this.makeContrast(probeAnalysisResult, factorsForName, term, factorName, contrastPvalue, interactionContrastTStats, interactionContrastCoeffs);
                    }
                } else {
                    if (!warned) {
                        LinearModelAnalyzer.log.warn("Interaction could not be computed for " + el + ", further warnings suppressed");
                        warned = true;
                    }
                    if (LinearModelAnalyzer.log.isDebugEnabled())
                        LinearModelAnalyzer.log.debug("Interaction could not be computed for " + el + ", further warnings suppressed");
                    // will over count?
                    notUsable++;
                    continue;
                }
            } else {
                /*
                     * Main effect
                     */
                assert factorsForName.size() == 1;
                ExperimentalFactor experimentalFactor = factorsForName.iterator().next();
                if (interceptFactor != null && factorsForName.size() == 1 && experimentalFactor.equals(interceptFactor)) {
                    overallPValue = lm.getInterceptP();
                } else {
                    overallPValue = lm.getMainEffectP(factorName);
                }
                /*
                     * Add contrasts unless overall pvalue is NaN
                     */
                if (overallPValue != null && !Double.isNaN(overallPValue)) {
                    Map<String, Double> mainEffectContrastTStats = lm.getContrastTStats(factorName);
                    Map<String, Double> mainEffectContrastPvalues = lm.getContrastPValues(factorName);
                    Map<String, Double> mainEffectContrastCoeffs = lm.getContrastCoefficients(factorName);
                    for (String term : mainEffectContrastPvalues.keySet()) {
                        Double contrastPvalue = mainEffectContrastPvalues.get(term);
                        this.makeContrast(probeAnalysisResult, factorsForName, term, factorName, contrastPvalue, mainEffectContrastTStats, mainEffectContrastCoeffs);
                    }
                } else {
                    if (!warned) {
                        LinearModelAnalyzer.log.warn("ANOVA could not be done for " + experimentalFactor + " on " + el + ", further warnings suppressed");
                        warned = true;
                    }
                    if (LinearModelAnalyzer.log.isDebugEnabled())
                        LinearModelAnalyzer.log.debug("ANOVA could not be done for " + experimentalFactor + " on " + el);
                    // will over count?
                    notUsable++;
                    continue;
                }
            }
            assert !Double.isNaN(overallPValue) : "We should not be keeping non-number pvalues (null or NaNs)";
            probeAnalysisResult.setPvalue(this.nan2Null(overallPValue));
            pvaluesForQvalue.get(factorName).add(overallPValue);
            resultLists.get(factorName).add(probeAnalysisResult);
        }
    // over terms
    }
    if (notUsable > 0) {
        LinearModelAnalyzer.log.info(notUsable + " elements or results were not usable - model could not be fit, etc.");
    }
    this.getRanksAndQvalues(resultLists, pvaluesForQvalue);
    DifferentialExpressionAnalysis expressionAnalysis = this.makeAnalysisEntity(bioAssaySet, config, label2Factors, baselineConditions, interceptFactor, interactionFactorLists, oneSampleTTest, resultLists, subsetFactorValue);
    LinearModelAnalyzer.log.info("Analysis processing phase done ...");
    return expressionAnalysis;
}