Examples with PI org.apache.commons.math3.util.FastMath.PI used on opensource projects

Example 6 with PI

use of org.apache.commons.math3.util.FastMath.PI in project knime-core by knime.

the class IrlsLearner method irlsRls.

/**
 * Do an irls step. The result is stored in beta.
 *
 * @param data over trainings data.
 * @param beta parameter vector
 * @param rC regressors count
 * @param tcC target category count
 * @throws CanceledExecutionException when method is cancelled
 */
private void irlsRls(final TrainingData<ClassificationTrainingRow> data, final RealMatrix beta, final int rC, final int tcC, final ExecutionMonitor exec) throws CanceledExecutionException {
    long rowCount = 0;
    int dim = (rC + 1) * (tcC - 1);
    RealMatrix xTwx = MatrixUtils.createRealMatrix(dim, dim);
    RealMatrix xTyu = MatrixUtils.createRealMatrix(dim, 1);
    double[] eBetaTx = new double[tcC - 1];
    double[] pi = new double[tcC - 1];
    final long totalRowCount = data.getRowCount();
    for (ClassificationTrainingRow row : data) {
        rowCount++;
        exec.checkCanceled();
        exec.setProgress(rowCount / (double) totalRowCount, "Row " + rowCount + "/" + totalRowCount);
        for (int k = 0; k < tcC - 1; k++) {
            double z = 0.0;
            for (FeatureIterator iter = row.getFeatureIterator(); iter.next(); ) {
                double featureVal = iter.getFeatureValue();
                int featureIdx = iter.getFeatureIndex();
                z += featureVal * beta.getEntry(0, k * (rC + 1) + featureIdx);
            }
            eBetaTx[k] = Math.exp(z);
        }
        double sumEBetaTx = 0;
        for (int k = 0; k < tcC - 1; k++) {
            sumEBetaTx += eBetaTx[k];
        }
        for (int k = 0; k < tcC - 1; k++) {
            double pik = eBetaTx[k] / (1 + sumEBetaTx);
            pi[k] = pik;
        }
        // fill xTwx (aka the hessian of the loglikelihood)
        for (FeatureIterator outer = row.getFeatureIterator(); outer.next(); ) {
            int i = outer.getFeatureIndex();
            double outerVal = outer.getFeatureValue();
            for (FeatureIterator inner = outer.spawn(); inner.next(); ) {
                int ii = inner.getFeatureIndex();
                double innerVal = inner.getFeatureValue();
                for (int k = 0; k < tcC - 1; k++) {
                    for (int kk = k; kk < tcC - 1; kk++) {
                        int o1 = k * (rC + 1);
                        int o2 = kk * (rC + 1);
                        double v = xTwx.getEntry(o1 + i, o2 + ii);
                        if (k == kk) {
                            double w = pi[k] * (1 - pi[k]);
                            v += outerVal * w * innerVal;
                            assert o1 == o2;
                        } else {
                            double w = -pi[k] * pi[kk];
                            v += outerVal * w * innerVal;
                        }
                        xTwx.setEntry(o1 + i, o2 + ii, v);
                        xTwx.setEntry(o1 + ii, o2 + i, v);
                        if (k != kk) {
                            xTwx.setEntry(o2 + ii, o1 + i, v);
                            xTwx.setEntry(o2 + i, o1 + ii, v);
                        }
                    }
                }
            }
        }
        int g = row.getCategory();
        // fill matrix xTyu
        for (FeatureIterator iter = row.getFeatureIterator(); iter.next(); ) {
            int idx = iter.getFeatureIndex();
            double val = iter.getFeatureValue();
            for (int k = 0; k < tcC - 1; k++) {
                int o = k * (rC + 1);
                double v = xTyu.getEntry(o + idx, 0);
                double y = k == g ? 1 : 0;
                v += (y - pi[k]) * val;
                xTyu.setEntry(o + idx, 0, v);
            }
        }
    }
    // currently not used but could become interesting in the future
    // if (m_penaltyTerm > 0.0) {
    // RealMatrix stdError = getStdErrorMatrix(xTwx);
    // // do not penalize the constant terms
    // for (int i = 0; i < tcC - 1; i++) {
    // stdError.setEntry(i * (rC + 1), i * (rC + 1), 0);
    // }
    // xTwx = xTwx.add(stdError.scalarMultiply(-0.00001));
    // }
    exec.checkCanceled();
    b = xTwx.multiply(beta.transpose()).add(xTyu);
    A = xTwx;
    if (rowCount < A.getColumnDimension()) {
        // but it's important to ensure this property
        throw new IllegalStateException("The dataset must have at least " + A.getColumnDimension() + " rows, but it has only " + rowCount + " rows. It is recommended to use a " + "larger dataset in order to increase accuracy.");
    }
    DecompositionSolver solver = new SingularValueDecomposition(A).getSolver();
    RealMatrix betaNew = solver.solve(b);
    beta.setSubMatrix(betaNew.transpose().getData(), 0, 0);
}

Example 7 with PI

use of org.apache.commons.math3.util.FastMath.PI in project ffx by mjschnie.

the class ReflectionList method getepsilon.

private void getepsilon(HKL hkl) {
    int epsilon = 1;
    int allowed = 255;
    int nsym = spaceGroup.symOps.size();
    for (int i = 1; i < nsym; i++) {
        HKL mate = new HKL();
        crystal.applySymRot(hkl, mate, spaceGroup.symOps.get(i));
        double shift = Crystal.sym_phase_shift(hkl, spaceGroup.symOps.get(i));
        if (mate.equals(hkl)) {
            if (cos(shift) > 0.999) {
                epsilon++;
            } else {
                allowed = 0;
                epsilon = 0;
                break;
            }
        } else if (mate.equals(HKL.neg(hkl))) {
            // centric reflection
            allowed = (int) rint(Crystal.mod(-0.5 * shift, PI) / (PI / HKL.ndiv));
        }
    }
    if (hkl.h() == 0 && hkl.k() == 0 && hkl.l() == 0) {
        allowed = 0;
    }
    hkl.epsilon(epsilon);
    hkl.allowed(allowed);
}

Example 8 with PI

use of org.apache.commons.math3.util.FastMath.PI in project tetrad by cmu-phil.

the class DataUtils method getNonparanormalTransformed.

// function (x, npn.func = "shrinkage", npn.thresh = NULL, verbose = TRUE)
// {
// gcinfo(FALSE)
// n = nrow(x)
// d = ncol(x)
// x.col = colnames(x)
// x.row = rownames(x)
// if (npn.func == "shrinkage") {
// if (verbose)
// cat("Conducting the nonparanormal (npn) transformation via shrunkun ECDF....")
// x = qnorm(apply(x, 2, rank)/(n + 1))
// x = x/sd(x[, 1])
// if (verbose)
// cat("done.\n")
// rm(n, d, verbose)
// gc()
// colnames(x) = x.col
// rownames(x) = x.row
// }
// if (npn.func == "truncation") {
// if (verbose)
// cat("Conducting nonparanormal (npn) transformation via truncated ECDF....")
// if (is.null(npn.thresh))
// npn.thresh = 1/(4 * (n^0.25) * sqrt(pi * log(n)))
// x = qnorm(pmin(pmax(apply(x, 2, rank)/n, npn.thresh),
// 1 - npn.thresh))
// x = x/sd(x[, 1])
// if (verbose)
// cat("done.\n")
// rm(n, d, npn.thresh, verbose)
// gc()
// colnames(x) = x.col
// rownames(x) = x.row
// }
// if (npn.func == "skeptic") {
// if (verbose)
// cat("Conducting nonparanormal (npn) transformation via skeptic....")
// x = 2 * sin(pi/6 * cor(x, method = "spearman"))
// if (verbose)
// cat("done.\n")
// rm(n, d, verbose)
// gc()
// colnames(x) = x.col
// rownames(x) = x.col
// }
// return(x)
// }
public static DataSet getNonparanormalTransformed(DataSet dataSet) {
    final TetradMatrix data = dataSet.getDoubleData();
    final TetradMatrix X = data.like();
    final double n = dataSet.getNumRows();
    final double delta = 1.0 / (4.0 * Math.pow(n, 0.25) * Math.sqrt(Math.PI * Math.log(n)));
    final NormalDistribution normalDistribution = new NormalDistribution();
    double std = Double.NaN;
    for (int j = 0; j < data.columns(); j++) {
        final double[] x1 = data.getColumn(j).toArray();
        double std1 = StatUtils.sd(x1);
        double mu1 = StatUtils.mean(x1);
        double[] x = ranks(data, x1);
        for (int i = 0; i < x.length; i++) {
            x[i] /= n;
            if (x[i] < delta)
                x[i] = delta;
            if (x[i] > (1. - delta))
                x[i] = 1. - delta;
            x[i] = normalDistribution.inverseCumulativeProbability(x[i]);
        }
        if (Double.isNaN(std)) {
            std = StatUtils.sd(x);
        }
        for (int i = 0; i < x.length; i++) {
            x[i] /= std;
            x[i] *= std1;
            x[i] += mu1;
        }
        X.assignColumn(j, new TetradVector(x));
    }
    return ColtDataSet.makeContinuousData(dataSet.getVariables(), X);
}