Example usage for java.lang Double POSITIVE_INFINITY

List of usage examples for java.lang Double POSITIVE_INFINITY

Introduction

In this page you can find the example usage for java.lang Double POSITIVE_INFINITY.

Prototype

double POSITIVE_INFINITY

To view the source code for java.lang Double POSITIVE_INFINITY.

Click Source Link

Document

A constant holding the positive infinity of type double .

Usage

From source file:com.opengamma.analytics.financial.model.finitedifference.PeacemanRachfordFiniteDifference2Db.java

@Override
public double[][] solve(final ConvectionDiffusion2DPDEDataBundle pdeData, final int tSteps, final int xSteps,
        final int ySteps, final double tMax, final BoundaryCondition2D xLowerBoundary,
        final BoundaryCondition2D xUpperBoundary, final BoundaryCondition2D yLowerBoundary,
        final BoundaryCondition2D yUpperBoundary, final Cube<Double, Double, Double, Double> freeBoundary) {

    final double dt = tMax / (tSteps);
    final double dx = (xUpperBoundary.getLevel() - xLowerBoundary.getLevel()) / (xSteps);
    final double dy = (yUpperBoundary.getLevel() - yLowerBoundary.getLevel()) / (ySteps);
    final double dtdx2 = dt / dx / dx;
    final double dtdx = dt / dx;
    final double dtdy2 = dt / dy / dy;
    final double dtdy = dt / dy;

    final double[][] v = new double[xSteps + 1][ySteps + 1];
    final double[][] vt = new double[xSteps + 1][ySteps + 1];

    final double[] x = new double[xSteps + 1];
    final double[] y = new double[ySteps + 1];

    final double[] q = new double[xSteps + 1];
    final double[] r = new double[ySteps + 1];
    final double[][] mx = new double[xSteps + 1][xSteps + 1];
    final double[][] my = new double[ySteps + 1][ySteps + 1];

    double currentX = 0;
    double currentY = 0;

    for (int j = 0; j <= ySteps; j++) {
        currentY = yLowerBoundary.getLevel() + j * dy;
        y[j] = currentY;//from  w w  w .  j  a  v a2 s . co  m
    }
    for (int i = 0; i <= xSteps; i++) {
        currentX = xLowerBoundary.getLevel() + i * dx;
        x[i] = currentX;
        for (int j = 0; j <= ySteps; j++) {
            v[i][j] = pdeData.getInitialValue(x[i], y[j]);
        }
    }

    double t, tStar;
    double a, b, c, d, f;

    for (int n = 0; n < tSteps; n++) {

        t = n * dt;
        tStar = t + 0.25 * dt;

        // stag 1 Explicit in y
        for (int i = 0; i <= xSteps; i++) {
            for (int j = 1; j < ySteps; j++) {
                c = pdeData.getC(t, x[i], y[j]);
                d = pdeData.getD(t, x[i], y[j]);
                f = pdeData.getF(t, x[i], y[j]);

                vt[i][j] = (1 - 0.25 * dt * c) * v[i][j];
                vt[i][j] -= 0.5 * dtdy2 * d * (v[i][j + 1] + v[i][j - 1] - 2 * v[i][j]);
                vt[i][j] -= 0.25 * dtdy * f * (v[i][j + 1] - v[i][j - 1]);
            }

            double[] temp = yLowerBoundary.getRightMatrixCondition(tStar, x[i]);
            double sum = 0;
            for (int k = 0; k < temp.length; k++) {
                sum += temp[k] * v[i][k];
            }
            sum += yLowerBoundary.getConstant(tStar, x[i], dy);

            temp = yLowerBoundary.getLeftMatrixCondition(tStar, x[i]);
            for (int k = 1; k < temp.length; k++) {
                sum -= temp[k] * vt[i][k];
            }
            vt[i][0] = sum / temp[0];

            temp = yUpperBoundary.getRightMatrixCondition(tStar, x[i]);
            sum = 0;
            for (int k = 0; k < temp.length; k++) {
                sum += temp[k] * v[i][ySteps - k];
            }
            sum += yUpperBoundary.getConstant(tStar, x[i], dy);

            temp = yUpperBoundary.getLeftMatrixCondition(tStar, x[i]);
            for (int k = 1; k < temp.length; k++) {
                sum -= temp[k] * vt[i][ySteps - k];
            }
            vt[i][ySteps] = sum / temp[0];
        }

        // stag 2 - Implicit in x
        t += 0.5 * dt;
        for (int j = 0; j <= ySteps; j++) {
            for (int i = 1; i < xSteps; i++) {
                a = pdeData.getA(t, x[i], y[j]);
                b = pdeData.getB(t, x[i], y[j]);
                c = pdeData.getC(t, x[i], y[j]);

                mx[i][i - 1] = 0.5 * (dtdx2 * a - 0.5 * dtdx * b);
                mx[i][i] = 1 + 0.5 * (-2 * dtdx2 * a + 0.5 * dt * c);
                mx[i][i + 1] = 0.5 * (dtdx2 * a + 0.5 * dtdx * b);

                q[i] = vt[i][j];
            }

            double[] temp = xLowerBoundary.getLeftMatrixCondition(t, y[j]);
            for (int k = 0; k < temp.length; k++) {
                mx[0][k] = temp[k];
            }
            temp = xUpperBoundary.getLeftMatrixCondition(t, y[j]);
            for (int k = 0; k < temp.length; k++) {
                mx[xSteps][xSteps - k] = temp[k];
            }

            temp = xLowerBoundary.getRightMatrixCondition(t, y[j]);
            double sum = 0;
            for (int k = 0; k < temp.length; k++) {
                sum += temp[k] * v[k][j];
            }
            q[0] = sum + xLowerBoundary.getConstant(t, y[j], dx);

            temp = xUpperBoundary.getRightMatrixCondition(t, y[j]);
            sum = 0;
            for (int k = 0; k < temp.length; k++) {
                sum += temp[k] * v[xSteps - k][j];
            }
            q[xSteps] = sum + xUpperBoundary.getConstant(t, y[j], dx);

            // SOR
            final double omega = 1.5;
            double scale = 1.0;
            double errorSqr = Double.POSITIVE_INFINITY;
            int min, max;
            int count = 0;
            while (errorSqr / (scale + 1e-10) > 1e-18 && count < SOR_MAX) {
                errorSqr = 0.0;
                scale = 0.0;
                for (int l = 0; l <= xSteps; l++) {
                    min = (l == xSteps ? 0 : Math.max(0, l - 1));
                    max = (l == 0 ? xSteps : Math.min(xSteps, l + 1));
                    sum = 0;
                    // for (int k = 0; k <= xSteps; k++) {
                    for (int k = min; k <= max; k++) { // mx is tri-diagonal so only need 3 steps here
                        sum += mx[l][k] * vt[k][j];
                    }
                    final double correction = omega / mx[l][l] * (q[l] - sum);
                    // if (freeBoundary != null) {
                    // correction = Math.max(correction, freeBoundary.getZValue(t, x[j]) - f[j]);
                    // }
                    errorSqr += correction * correction;
                    vt[l][j] += correction;
                    scale += vt[l][j] * vt[l][j];
                }
                count++;
            }
            Validate.isTrue(count < SOR_MAX, "SOR exceeded max interations");
        }

        // stag 3 explicit in x
        tStar = t + 0.25 * dt;
        for (int j = 0; j <= ySteps; j++) {
            for (int i = 1; i < xSteps; i++) {

                a = pdeData.getA(t, x[i], y[j]);
                b = pdeData.getB(t, x[i], y[j]);
                c = pdeData.getC(t, x[i], y[j]);

                v[i][j] = (1 - 0.25 * c) * vt[i][j];
                v[i][j] -= 0.5 * dtdx2 * a * (vt[i + 1][j] + vt[i - 1][j] - 2 * vt[i][j]);
                v[i][j] -= 0.25 * dtdx * b * (vt[i + 1][j] - vt[i - 1][j]);
            }

            double[] temp = xLowerBoundary.getRightMatrixCondition(tStar, y[j]);
            double sum = 0;
            for (int k = 0; k < temp.length; k++) {
                sum += temp[k] * vt[k][j];
            }
            sum += xLowerBoundary.getConstant(tStar, y[j], dx);

            temp = xLowerBoundary.getLeftMatrixCondition(tStar, y[j]);
            for (int k = 1; k < temp.length; k++) {
                sum -= temp[k] * v[k][j];
            }
            v[0][j] = sum / temp[0];

            temp = xUpperBoundary.getRightMatrixCondition(tStar, y[j]);
            sum = 0;
            for (int k = 0; k < temp.length; k++) {
                sum += temp[k] * vt[xSteps - k][j];
            }
            sum += xUpperBoundary.getConstant(tStar, y[j], dx);

            temp = xUpperBoundary.getLeftMatrixCondition(tStar, y[j]);
            for (int k = 1; k < temp.length; k++) {
                sum -= temp[k] * v[xSteps - k][j];
            }
            v[xSteps][j] = sum / temp[0];
        }

        // stag 4 - implicit in y
        t = (n + 1) * dt;
        for (int i = 0; i <= xSteps; i++) {
            for (int j = 1; j < ySteps; j++) {

                c = pdeData.getC(t, x[i], y[j]);
                d = pdeData.getD(t, x[i], y[j]);
                f = pdeData.getF(t, x[i], y[j]);

                my[j][j - 1] = 0.5 * (dtdy2 * d - 0.5 * dtdy * f);
                my[j][j] = 1 + 0.5 * (-2 * dtdy2 * d + 0.5 * dt * c);
                my[j][j + 1] = 0.5 * (dtdy2 * d + 0.5 * dtdy * f);

                r[j] = v[i][j];
            }

            double[] temp = yLowerBoundary.getLeftMatrixCondition(t, x[i]);
            for (int k = 0; k < temp.length; k++) {
                my[0][k] = temp[k];
            }
            temp = yUpperBoundary.getLeftMatrixCondition(t, x[i]);
            for (int k = 0; k < temp.length; k++) {
                my[ySteps][ySteps - k] = temp[k];
            }

            temp = yLowerBoundary.getRightMatrixCondition(t, x[i]);
            double sum = 0;
            for (int k = 0; k < temp.length; k++) {
                sum += temp[k] * vt[i][k];
            }
            r[0] = sum + yLowerBoundary.getConstant(t, x[i], dy);

            temp = yUpperBoundary.getRightMatrixCondition(t, x[i]);
            sum = 0;
            for (int k = 0; k < temp.length; k++) {
                sum += temp[k] * vt[i][ySteps - k];
            }
            r[ySteps] = sum + yUpperBoundary.getConstant(t, x[i], dy);

            // SOR
            final double omega = 1.5;
            double scale = 1.0;
            double errorSqr = Double.POSITIVE_INFINITY;
            int count = 0;
            while (errorSqr / (scale + 1e-10) > 1e-18 && count < SOR_MAX) {
                errorSqr = 0.0;
                scale = 0.0;
                int min, max;
                for (int l = 0; l <= ySteps; l++) {
                    min = (l == ySteps ? 0 : Math.max(0, l - 1));
                    max = (l == 0 ? ySteps : Math.min(ySteps, l + 1));
                    sum = 0;
                    // for (int k = 0; k <= ySteps; k++) {
                    for (int k = min; k <= max; k++) {
                        sum += my[l][k] * v[i][k];
                    }
                    final double correction = omega / my[l][l] * (r[l] - sum);
                    // if (freeBoundary != null) {
                    // correction = Math.max(correction, freeBoundary.getZValue(t, x[j]) - f[j]);
                    // }
                    errorSqr += correction * correction;
                    v[i][l] += correction;
                    scale += v[i][l] * v[i][l];
                }
                count++;
            }
            Validate.isTrue(count < SOR_MAX, "SOR exceeded max interations");
        }

    } // time loop
    return v;

}

From source file:iDynoOptimizer.MOEAFramework26.src.org.moeaframework.analysis.diagnostics.ApproximationSetViewer.java

/**
 * Initializes the reference set bounds.
 *//*from w  ww .j a v a 2s  .  com*/
protected void initializeReferenceSetBounds() {
    if (referenceSet == null) {
        return;
    }

    double domainMin = Double.POSITIVE_INFINITY;
    double domainMax = Double.NEGATIVE_INFINITY;
    double rangeMin = Double.POSITIVE_INFINITY;
    double rangeMax = Double.NEGATIVE_INFINITY;

    for (Solution solution : referenceSet) {
        domainMin = Math.min(domainMin, getValue(solution, 0));
        domainMax = Math.max(domainMax, getValue(solution, 0));
        rangeMin = Math.min(rangeMin, getValue(solution, 1));
        rangeMax = Math.max(rangeMax, getValue(solution, 1));
    }

    domainMax += (domainMax - domainMin);
    rangeMax += (rangeMax - rangeMin);

    referenceDomainBounds = new Range(domainMin, domainMax);
    referenceRangeBounds = new Range(rangeMin, rangeMax);
}

From source file:it.unibo.alchemist.model.implementations.reactions.SAPEREReaction.java

@Override
public void execute() {
    if (possibleMatches.isEmpty()) {
        for (final ILsaAction a : getSAPEREActions()) {
            a.setExecutionContext(null, validNodes);
            a.execute();/*from www.j  a va  2 s . c o  m*/
        }
        return;
    }
    final Position nodePosCache = modifiesOnlyLocally ? environment.getPosition(getNode()) : null;
    final List<? extends ILsaMolecule> localContentCache = modifiesOnlyLocally
            ? new ArrayList<>(getNode().getLsaSpace())
            : null;
    Map<HashString, ITreeNode<?>> matches = null;
    Map<ILsaNode, List<ILsaMolecule>> toRemove = null;
    /*
     * If there is infinite propensity, the last match added is the one to
     * choose, since it is the one which generated the "infinity" value.
     */
    if (totalPropensity == Double.POSITIVE_INFINITY) {
        final int index = possibleMatches.size() - 1;
        matches = possibleMatches.get(index);
        toRemove = possibleRemove.get(index);
    } else if (numericRate()) {
        /*
         * If the rate is numeric, the choice is just random
         */
        final int index = Math.abs(rng.nextInt()) % possibleMatches.size();
        matches = possibleMatches.get(index);
        toRemove = possibleRemove.get(index);
    } else {
        /*
         * Otherwise, the matches must be chosen randomly using their
         * propensities
         */
        final double index = rng.nextDouble() * totalPropensity;
        double sum = 0;
        for (int i = 0; matches == null; i++) {
            sum += propensities.get(i);
            if (sum > index) {
                matches = possibleMatches.get(i);
                toRemove = possibleRemove.get(i);
            }
        }
    }
    /*
     * The matched LSAs must be removed from the local space, if no action
     * added them back.
     */
    for (final Entry<ILsaNode, List<ILsaMolecule>> entry : toRemove.entrySet()) {
        final ILsaNode n = entry.getKey();
        for (final ILsaMolecule m : entry.getValue()) {
            n.removeConcentration(m);
        }
    }
    /*
     * #T Must be loaded by the reaction, which is the only structure aware
     * of the time. Other special values (#NEIG, #O, #D) will be allocated
     * inside the actions.
     */
    matches.put(LsaMolecule.SYN_T, new NumTreeNode(getTau().toDouble()));
    for (final ILsaAction a : getSAPEREActions()) {
        a.setExecutionContext(matches, validNodes);
        a.execute();
    }

    /*
     * Empty action optimization
     */
    if (modifiesOnlyLocally) {
        final ILsaNode n = getNode();
        if (nodePosCache.equals(environment.getPosition(getNode()))) {
            final List<? extends ILsaMolecule> contents = n.getLsaSpace();
            if (contents.size() == localContentCache.size()) {
                emptyExecution = localContentCache.containsAll(contents);
            }
        }
    }
}

From source file:com.rapidminer.operator.preprocessing.discretization.AbsoluteDiscretization.java

@Override
public PreprocessingModel createPreprocessingModel(ExampleSet exampleSet) throws OperatorException {
    DiscretizationModel model = new DiscretizationModel(exampleSet);

    exampleSet.recalculateAllAttributeStatistics();
    // calculating number of bins
    int sizeOfBins = getParameterAsInt(PARAMETER_SIZE_OF_BINS);
    int numberOfBins = exampleSet.size() / sizeOfBins;
    int numberOfExamples = exampleSet.size();
    // add one bin if a remainder exists
    if (numberOfBins * sizeOfBins < numberOfExamples) {
        numberOfBins++;/*  w ww .j a va  2 s .  c o m*/
    }
    HashMap<Attribute, double[]> ranges = new HashMap<Attribute, double[]>();

    int sortingDirection = getParameterAsInt(PARAMETER_SORTING_DIRECTION);

    for (Attribute attribute : exampleSet.getAttributes()) {
        if (attribute.isNumerical()) { // skip nominal and date attributes
            ExampleSet sortedSet = new SortedExampleSet(exampleSet, attribute, sortingDirection);
            double[] binRange = new double[numberOfBins];
            for (int i = 0; i < numberOfBins - 1; i++) {
                int offset = (i + 1) * sizeOfBins - 1;
                double infimum = sortedSet.getExample(offset).getValue(attribute);
                offset++;
                double supremum = sortedSet.getExample(offset).getValue(attribute);
                // if targets equal values: Search for next different value
                while (infimum == supremum && offset < numberOfExamples) {
                    supremum = sortedSet.getExample(offset).getValue(attribute);
                    offset++;
                }
                if (sortingDirection == SortedExampleSet.DECREASING) {
                    binRange[numberOfBins - 2 - i] = (infimum + supremum) / 2d;
                } else {
                    binRange[i] = (infimum + supremum) / 2d;
                }
            }
            binRange[numberOfBins - 1] = Double.POSITIVE_INFINITY;
            ranges.put(attribute, binRange);
        }
    }

    // determine number of digits
    int numberOfDigits = -1;
    if (getParameterAsBoolean(PARAMETER_AUTOMATIC_NUMBER_OF_DIGITS) == false) {
        numberOfDigits = getParameterAsInt(PARAMETER_NUMBER_OF_DIGITS);
    }

    model.setRanges(ranges, "range", getParameterAsInt(PARAMETER_RANGE_NAME_TYPE), numberOfDigits);
    return (model);
}

From source file:maltcms.ui.fileHandles.csv.CSV2JFCLoader.java

@Override
public void run() {
    CSV2TableLoader tl = new CSV2TableLoader(this.ph, this.is);

    DefaultTableModel dtm;// www. j  a va 2 s  .  com
    try {
        dtm = tl.call();
        if (this.mode == CHART.XY) {
            XYSeriesCollection cd = new XYSeriesCollection();
            for (int j = 0; j < dtm.getColumnCount(); j++) {
                XYSeries xys = new XYSeries(dtm.getColumnName(j));
                for (int i = 0; i < dtm.getRowCount(); i++) {
                    Object o = dtm.getValueAt(i, j);
                    try {
                        double d = Double.parseDouble(o.toString());
                        xys.add(i, d);
                        Logger.getLogger(getClass().getName()).log(Level.INFO, "Adding {0} {1} {2}",
                                new Object[] { i, d, dtm.getColumnName(j) });
                    } catch (Exception e) {
                    }
                }
                cd.addSeries(xys);
            }
            XYLineAndShapeRenderer d = new XYLineAndShapeRenderer(true, false);
            XYPlot xyp = new XYPlot(cd, new NumberAxis("category"), new NumberAxis("value"), d);

            JFreeChart jfc = new JFreeChart(this.title, xyp);
            jtc.setChart(jfc);
            Logger.getLogger(getClass().getName()).info("creating chart done");
        } else if (this.mode == CHART.MATRIX) {
            DefaultXYZDataset cd = new DefaultXYZDataset();
            Logger.getLogger(getClass().getName()).log(Level.INFO, "Name of column 0: {0}",
                    dtm.getColumnName(0));
            if (dtm.getColumnName(0).isEmpty()) {
                Logger.getLogger(getClass().getName()).info("Removing column 0");
                dtm = removeColumn(dtm, 0);
            }
            if (dtm.getColumnName(dtm.getColumnCount() - 1).equalsIgnoreCase("filename")) {
                dtm = removeColumn(dtm, dtm.getColumnCount() - 1);
            }
            StringBuilder sb = new StringBuilder();
            for (int i = 0; i < dtm.getRowCount(); i++) {
                for (int j = 0; j < dtm.getColumnCount(); j++) {
                    sb.append(dtm.getValueAt(i, j) + " ");
                }
                sb.append("\n");
            }
            Logger.getLogger(getClass().getName()).log(Level.INFO, "Table before sorting: {0}", sb.toString());
            //                dtm = sort(dtm);
            StringBuilder sb2 = new StringBuilder();
            for (int i = 0; i < dtm.getRowCount(); i++) {
                for (int j = 0; j < dtm.getColumnCount(); j++) {
                    sb2.append(dtm.getValueAt(i, j) + " ");
                }
                sb2.append("\n");
            }
            Logger.getLogger(getClass().getName()).log(Level.INFO, "Table after sorting: {0}", sb2.toString());
            int rows = dtm.getRowCount();
            int columns = dtm.getColumnCount();
            Logger.getLogger(getClass().getName()).log(Level.INFO, "Storing {0} * {1} elements, {2} total!",
                    new Object[] { columns, rows, rows * columns });
            double[][] data = new double[3][(columns * rows)];
            ArrayDouble.D1 dt = new ArrayDouble.D1((columns) * rows);
            double min = Double.POSITIVE_INFINITY;
            double max = Double.NEGATIVE_INFINITY;
            EvalTools.eqI(rows, columns, this);
            int k = 0;
            for (int i = 0; i < dtm.getRowCount(); i++) {
                for (int j = 0; j < dtm.getColumnCount(); j++) {

                    Object o = dtm.getValueAt(i, j);
                    try {
                        double d = Double.parseDouble(o.toString());
                        if (d < min) {
                            min = d;
                        }
                        if (d > max) {
                            max = d;
                        }
                        data[0][k] = (double) i;
                        data[1][k] = (double) j;
                        data[2][k] = d;
                        dt.set(k, d);
                        k++;
                        //System.out.println("Adding "+i+" "+d+" "+dtm.getColumnName(j));
                    } catch (Exception e) {
                    }
                }
                //cd.addSeries(xys);
            }
            cd.addSeries(this.title, data);
            XYBlockRenderer xyb = new XYBlockRenderer();
            GradientPaintScale ps = new GradientPaintScale(ImageTools.createSampleTable(256), min, max,
                    ImageTools
                            .rampToColorArray(new ColorRampReader().readColorRamp("res/colorRamps/bcgyr.csv")));

            xyb.setPaintScale(ps);
            final String[] colnames = new String[dtm.getColumnCount()];
            for (int i = 0; i < colnames.length; i++) {
                colnames[i] = dtm.getColumnName(i);
            }
            NumberAxis na1 = new SymbolAxis("category", colnames);
            na1.setVerticalTickLabels(false);
            NumberAxis na2 = new SymbolAxis("category", colnames);
            na1.setVerticalTickLabels(true);
            XYPlot xyp = new XYPlot(cd, na1, na2, xyb);
            xyb.setSeriesToolTipGenerator(0, new XYToolTipGenerator() {

                @Override
                public String generateToolTip(XYDataset xyd, int i, int i1) {
                    return "[" + colnames[xyd.getX(i, i1).intValue()] + ":"
                            + colnames[xyd.getY(i, i1).intValue()] + "] = "
                            + ((XYZDataset) xyd).getZValue(i, i1) + "";
                }
            });

            JFreeChart jfc = new JFreeChart(this.title, xyp);
            NumberAxis values = new NumberAxis("value");
            values.setAutoRange(false);
            values.setRangeWithMargins(min, max);
            PaintScaleLegend psl = new PaintScaleLegend(ps, values);
            psl.setBackgroundPaint(jfc.getBackgroundPaint());
            jfc.addSubtitle(psl);
            psl.setStripWidth(50);
            psl.setPadding(20, 20, 20, 20);
            psl.setHeight(200);
            psl.setPosition(RectangleEdge.RIGHT);
            jtc.setChart(jfc);
        }
    } catch (Exception ex) {
        Exceptions.printStackTrace(ex);
    }

    ph.finish();
}

From source file:hivemall.utils.math.MathUtils.java

/**
 * Returns the inverse erf. This code is based on erfInv() in
 * org.apache.commons.math3.special.Erf.
 * <p>/*  w  w  w. j  a  va2 s.com*/
 * This implementation is described in the paper: <a
 * href="http://people.maths.ox.ac.uk/gilesm/files/gems_erfinv.pdf">Approximating the erfinv
 * function</a> by Mike Giles, Oxford-Man Institute of Quantitative Finance, which was published
 * in GPU Computing Gems, volume 2, 2010. The source code is available <a
 * href="http://gpucomputing.net/?q=node/1828">here</a>.
 * </p>
 * 
 * @param x the value
 * @return t such that x = erf(t)
 */
public static double inverseErf(final double x) {

    // beware that the logarithm argument must be
    // commputed as (1.0 - x) * (1.0 + x),
    // it must NOT be simplified as 1.0 - x * x as this
    // would induce rounding errors near the boundaries +/-1
    double w = -Math.log((1.0 - x) * (1.0 + x));
    double p;

    if (w < 6.25) {
        w = w - 3.125;
        p = -3.6444120640178196996e-21;
        p = -1.685059138182016589e-19 + p * w;
        p = 1.2858480715256400167e-18 + p * w;
        p = 1.115787767802518096e-17 + p * w;
        p = -1.333171662854620906e-16 + p * w;
        p = 2.0972767875968561637e-17 + p * w;
        p = 6.6376381343583238325e-15 + p * w;
        p = -4.0545662729752068639e-14 + p * w;
        p = -8.1519341976054721522e-14 + p * w;
        p = 2.6335093153082322977e-12 + p * w;
        p = -1.2975133253453532498e-11 + p * w;
        p = -5.4154120542946279317e-11 + p * w;
        p = 1.051212273321532285e-09 + p * w;
        p = -4.1126339803469836976e-09 + p * w;
        p = -2.9070369957882005086e-08 + p * w;
        p = 4.2347877827932403518e-07 + p * w;
        p = -1.3654692000834678645e-06 + p * w;
        p = -1.3882523362786468719e-05 + p * w;
        p = 0.0001867342080340571352 + p * w;
        p = -0.00074070253416626697512 + p * w;
        p = -0.0060336708714301490533 + p * w;
        p = 0.24015818242558961693 + p * w;
        p = 1.6536545626831027356 + p * w;
    } else if (w < 16.0) {
        w = Math.sqrt(w) - 3.25;
        p = 2.2137376921775787049e-09;
        p = 9.0756561938885390979e-08 + p * w;
        p = -2.7517406297064545428e-07 + p * w;
        p = 1.8239629214389227755e-08 + p * w;
        p = 1.5027403968909827627e-06 + p * w;
        p = -4.013867526981545969e-06 + p * w;
        p = 2.9234449089955446044e-06 + p * w;
        p = 1.2475304481671778723e-05 + p * w;
        p = -4.7318229009055733981e-05 + p * w;
        p = 6.8284851459573175448e-05 + p * w;
        p = 2.4031110387097893999e-05 + p * w;
        p = -0.0003550375203628474796 + p * w;
        p = 0.00095328937973738049703 + p * w;
        p = -0.0016882755560235047313 + p * w;
        p = 0.0024914420961078508066 + p * w;
        p = -0.0037512085075692412107 + p * w;
        p = 0.005370914553590063617 + p * w;
        p = 1.0052589676941592334 + p * w;
        p = 3.0838856104922207635 + p * w;
    } else if (!Double.isInfinite(w)) {
        w = Math.sqrt(w) - 5.0;
        p = -2.7109920616438573243e-11;
        p = -2.5556418169965252055e-10 + p * w;
        p = 1.5076572693500548083e-09 + p * w;
        p = -3.7894654401267369937e-09 + p * w;
        p = 7.6157012080783393804e-09 + p * w;
        p = -1.4960026627149240478e-08 + p * w;
        p = 2.9147953450901080826e-08 + p * w;
        p = -6.7711997758452339498e-08 + p * w;
        p = 2.2900482228026654717e-07 + p * w;
        p = -9.9298272942317002539e-07 + p * w;
        p = 4.5260625972231537039e-06 + p * w;
        p = -1.9681778105531670567e-05 + p * w;
        p = 7.5995277030017761139e-05 + p * w;
        p = -0.00021503011930044477347 + p * w;
        p = -0.00013871931833623122026 + p * w;
        p = 1.0103004648645343977 + p * w;
        p = 4.8499064014085844221 + p * w;
    } else {
        // this branch does not appears in the original code, it
        // was added because the previous branch does not handle
        // x = +/-1 correctly. In this case, w is positive infinity
        // and as the first coefficient (-2.71e-11) is negative.
        // Once the first multiplication is done, p becomes negative
        // infinity and remains so throughout the polynomial evaluation.
        // So the branch above incorrectly returns negative infinity
        // instead of the correct positive infinity.
        p = Double.POSITIVE_INFINITY;
    }

    return p * x;
}

From source file:edu.cmu.tetrad.util.StatUtils.java

/**
 * @param array a double array.//from w ww.j  a v  a 2s .c om
 * @param N     the number of values of array which should be considered.
 * @return the median of the first N values in this array.
 */
public static double median(double array[], int N) {

    double a[] = new double[N + 1];

    System.arraycopy(array, 0, a, 0, N);

    a[N] = Double.POSITIVE_INFINITY;

    double v, t;
    int i, j, l = 0;
    int r = N - 1, k1 = r / 2, k2 = r - k1;

    while (r > l) {
        v = a[l];
        i = l;
        j = r + 1;

        for (;;) {
            while (a[++i] < v) {
            }
            while (a[--j] > v) {
            }

            if (i >= j) {
                break;
            }

            t = a[i];
            a[i] = a[j];
            a[j] = t;
        }

        t = a[j];
        a[j] = a[l];
        a[l] = t;

        if (j <= k1) {
            l = j + 1;
        }

        if (j >= k2) {
            r = j - 1;
        }
    }

    return (a[k1] + a[k2]) / 2;
}

From source file:com.opengamma.analytics.financial.model.finitedifference.CrankNicolsonFiniteDifferenceSOR.java

@Override
public PDEResults1D solve(ConvectionDiffusionPDE1DStandardCoefficients pdeData,
        Function1D<Double, Double> initialCondition, int tSteps, int xSteps, double tMax,
        BoundaryCondition lowerBoundary, BoundaryCondition upperBoundary,
        Surface<Double, Double, Double> freeBoundary) {
    // simple test code - doesn't use a PDEGrid1D
    final PDEGrid1D grid = new PDEGrid1D(tSteps + 1, xSteps + 1, tMax, lowerBoundary.getLevel(),
            upperBoundary.getLevel());/*from   w w  w .  j a v  a 2 s. c  o m*/

    final double dt = tMax / (tSteps);
    final double dx = (upperBoundary.getLevel() - lowerBoundary.getLevel()) / (xSteps);
    final double dtdx2 = dt / dx / dx;
    final double dtdx = dt / dx;

    final double[] f = new double[xSteps + 1];
    final double[] x = new double[xSteps + 1];
    final double[] q = new double[xSteps + 1];
    final double[][] m = new double[xSteps + 1][xSteps + 1];
    // double[] coefficients = new double[3];

    double currentX = lowerBoundary.getLevel();

    double a, b, c, aa, bb, cc;

    for (int i = 0; i <= xSteps; i++) {
        currentX = lowerBoundary.getLevel() + i * dx;
        x[i] = currentX;
        final double value = initialCondition.evaluate(currentX);
        f[i] = value;
    }

    double t = 0.0;

    for (int n = 0; n < tSteps; n++) {
        t += dt;

        for (int i = 1; i < xSteps; i++) {
            a = pdeData.getA(t - dt, x[i]);
            b = pdeData.getB(t - dt, x[i]);
            c = pdeData.getC(t - dt, x[i]);

            double rho = a;
            // double bdx = (b * dx / 2);
            // if (Math.abs(bdx) > 10 * Math.abs(a)) {
            // rho = Math.abs(bdx);
            // } else if (Math.abs(a) > 10 * Math.abs(bdx)) {
            // rho = a;
            // } else {
            // rho = bdx / Math.tanh(bdx / a);
            // }

            aa = (1 - _theta) * (-dtdx2 * rho + 0.5 * dtdx * b);
            bb = 1 + (1 - _theta) * (2 * dtdx2 * rho - dt * c);
            cc = (1 - _theta) * (-dtdx2 * rho - 0.5 * dtdx * b);
            q[i] = aa * f[i - 1] + bb * f[i] + cc * f[i + 1];

            // TODO could store these
            a = pdeData.getA(t, x[i]);
            b = pdeData.getB(t, x[i]);
            c = pdeData.getC(t, x[i]);

            rho = a;
            // bdx = (b * dx / 2);
            // if (Math.abs(bdx) > 10 * Math.abs(a)) {
            // rho = Math.abs(bdx);
            // } else if (Math.abs(a) > 10 * Math.abs(bdx)) {
            // rho = a;
            // } else {
            // rho = bdx / Math.tanh(bdx / a);
            // }

            aa = (-dtdx2 * rho + 0.5 * dtdx * b);
            bb = (2 * dtdx2 * rho - dt * c);
            cc = (-dtdx2 * rho - 0.5 * dtdx * b);
            m[i][i - 1] = -_theta * aa;
            m[i][i] = 1 - _theta * bb;
            m[i][i + 1] = -_theta * cc;
        }

        double[] temp = lowerBoundary.getLeftMatrixCondition(pdeData, grid, t);
        for (int k = 0; k < temp.length; k++) {
            m[0][k] = temp[k];
        }
        temp = upperBoundary.getLeftMatrixCondition(pdeData, grid, t);
        for (int k = 0; k < temp.length; k++) {
            m[xSteps][xSteps - k] = temp[k];
        }

        temp = lowerBoundary.getRightMatrixCondition(pdeData, grid, t);
        double sum = 0;
        for (int k = 0; k < temp.length; k++) {
            sum += temp[k] * f[k];
        }
        q[0] = sum + lowerBoundary.getConstant(pdeData, t);

        temp = upperBoundary.getRightMatrixCondition(pdeData, grid, t);
        sum = 0;
        for (int k = 0; k < temp.length; k++) {
            sum += temp[k] * f[xSteps - k];
        }
        q[xSteps] = sum + upperBoundary.getConstant(pdeData, t);

        // SOR
        final double omega = 1.0;
        double scale = 1.0;
        double errorSqr = Double.POSITIVE_INFINITY;
        while (errorSqr / (scale + 1e-10) > 1e-18) {
            errorSqr = 0.0;
            scale = 0.0;
            for (int j = 0; j <= xSteps; j++) {
                sum = 0;
                for (int k = 0; k <= xSteps; k++) {
                    sum += m[j][k] * f[k];
                }
                double correction = omega / m[j][j] * (q[j] - sum);
                if (freeBoundary != null) {
                    correction = Math.max(correction, freeBoundary.getZValue(t, x[j]) - f[j]);
                }
                errorSqr += correction * correction;
                f[j] += correction;
                scale += f[j] * f[j];
            }
        }

    }

    return new PDETerminalResults1D(grid, f);
}

From source file:gdsc.smlm.fitting.BinomialFitter.java

/**
 * Fit the binomial distribution (n,p) to the input data. Performs fitting assuming a fixed n value and attempts to
 * optimise p. All n from minN to maxN are evaluated. If maxN is zero then all possible n from minN are evaluated
 * until the fit is worse.//  www. j  a v a  2 s . c om
 * 
 * @param data
 *            The input data (all value must be positive)
 * @param minN
 *            The minimum n to evaluate
 * @param maxN
 *            The maximum n to evaluate. Set to zero to evaluate all possible values.
 * @param zeroTruncated
 *            True if the model should ignore n=0 (zero-truncated binomial)
 * @return The best fit (n, p)
 * @throws IllegalArgumentException
 *             If any of the input data values are negative
 */
public double[] fitBinomial(int[] data, int minN, int maxN, boolean zeroTruncated) {
    double[] histogram = getHistogram(data, false);

    final double initialSS = Double.POSITIVE_INFINITY;
    double bestSS = initialSS;
    double[] parameters = null;
    int worse = 0;
    int N = (int) histogram.length - 1;
    if (minN < 1)
        minN = 1;
    if (maxN > 0) {
        if (N > maxN) {
            // Limit the number fitted to maximum
            N = maxN;
        } else if (N < maxN) {
            // Expand the histogram to the maximum
            histogram = Arrays.copyOf(histogram, maxN + 1);
            N = maxN;
        }
    }
    if (minN > N)
        minN = N;

    final double mean = getMean(histogram);

    String name = (zeroTruncated) ? "Zero-truncated Binomial distribution" : "Binomial distribution";

    log("Mean cluster size = %s", Utils.rounded(mean));
    log("Fitting cumulative " + name);

    // Since varying the N should be done in integer steps do this
    // for n=1,2,3,... until the SS peaks then falls off (is worse than the best 
    // score several times in succession)
    for (int n = minN; n <= N; n++) {
        PointValuePair solution = fitBinomial(histogram, mean, n, zeroTruncated);
        if (solution == null)
            continue;

        double p = solution.getPointRef()[0];

        log("Fitted %s : N=%d, p=%s. SS=%g", name, n, Utils.rounded(p), solution.getValue());

        if (bestSS > solution.getValue()) {
            bestSS = solution.getValue();
            parameters = new double[] { n, p };
            worse = 0;
        } else if (bestSS != initialSS) {
            if (++worse >= 3)
                break;
        }
    }

    return parameters;
}

From source file:hivemall.knn.similarity.DIMSUMMapperUDTF.java

@SuppressWarnings("unchecked")
@Override/*from   ww w .j a  va  2s  .c  om*/
public void process(Object[] args) throws HiveException {
    Feature[] row = parseFeatures(args[0]);
    if (row == null) {
        return;
    }
    this.probes = row;

    // since the 2nd argument (column norms) is consistent,
    // column-related values, `colNorms` and `colProbs`, should be cached
    if (colNorms == null || colProbs == null) {
        final int numCols = colNormsOI.getMapSize(args[1]);

        if (sqrtGamma == Double.POSITIVE_INFINITY) { // set default value to `gamma` based on `threshold`
            if (threshold > 0.d) { // if `threshold` = 0, `gamma` is INFINITY i.e. always accept <j, k> pairs
                this.sqrtGamma = Math.sqrt(10 * Math.log(numCols) / threshold);
            }
        }

        this.colNorms = new HashMap<Object, Double>(numCols);
        this.colProbs = new HashMap<Object, Double>(numCols);
        final Map<Object, Object> m = (Map<Object, Object>) colNormsOI.getMap(args[1]);
        for (Map.Entry<Object, Object> e : m.entrySet()) {
            Object j = e.getKey();
            if (parseFeatureAsInt) {
                j = HiveUtils.asJavaInt(j);
            } else {
                j = j.toString();
            }

            double norm = HiveUtils.asJavaDouble(e.getValue());
            if (norm == 0.d) { // avoid zero-division
                norm = 1.d;
            }

            colNorms.put(j, norm);

            double p = Math.min(1.d, sqrtGamma / norm);
            colProbs.put(j, p);
        }
    }

    if (parseFeatureAsInt) {
        forwardAsIntFeature(row);
    } else {
        forwardAsStringFeature(row);
    }
}