List of usage examples for java.lang Double POSITIVE_INFINITY
double POSITIVE_INFINITY
To view the source code for java.lang Double POSITIVE_INFINITY.
Click Source Link
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); } }