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:jp.furplag.util.commons.NumberUtilsTest.java
/** * {@link jp.furplag.util.commons.NumberUtils#add(java.lang.Number, java.lang.Number)}. *//* w w w.j a v a 2 s. c o m*/ @SuppressWarnings("unchecked") @Test public void testAddTNumber() { assertEquals("null", null, add(null, null)); assertEquals("null", null, add(null, 10)); assertEquals("null", (Object) 123.456d, add(123.456, 0)); assertEquals("null", (Object) 123.456f, add(123.456f, 0)); assertEquals("123 + .456: Float", (Object) 123, add(123, .456f)); assertEquals("123 + .456: Float", (Object) Double.class, add(123d, .456f).getClass()); assertEquals("123 + .456: Float", (Object) BigDecimal.class, add(BigDecimal.valueOf(123d), .456f).getClass()); for (Class<?> type : NUMBERS) { try { Object o = null; Number augend = null; Class<? extends Number> typeOfN = (Class<? extends Number>) type; Class<? extends Number> wrapper = (Class<? extends Number>) ClassUtils.primitiveToWrapper(type); if (ClassUtils.isPrimitiveWrapper(wrapper)) { o = wrapper.getMethod("valueOf", String.class).invoke(null, "123"); } else { Constructor<?> c = type.getDeclaredConstructor(String.class); o = c.newInstance("123"); } assertEquals("123.456: " + type.getSimpleName(), add((Number) o, augend, typeOfN), add((Number) o, augend)); augend = .456f; assertEquals("123.456: " + type.getSimpleName(), add((Number) o, augend, typeOfN), add((Number) o, augend)); augend = Double.NaN; assertEquals("NaN: " + type.getSimpleName(), add((Number) o, augend, typeOfN), add((Number) o, augend)); augend = Float.NEGATIVE_INFINITY; assertEquals("-Infinity: Float: " + type.getSimpleName(), add((Number) o, augend, typeOfN), add((Number) o, augend)); augend = Double.POSITIVE_INFINITY; assertEquals("Infinity: Double: : " + type.getSimpleName(), add((Number) o, augend, typeOfN), add((Number) o, augend)); } catch (Exception e) { e.printStackTrace(); fail(e.getMessage() + "\n" + Arrays.toString(e.getStackTrace())); } } }
From source file:ddf.catalog.source.opensearch.impl.OpenSearchParserImpl.java
/** * Takes in an array of coordinates and converts it to a (rough approximation) bounding box. * * <p>Note: Searches being performed where the polygon goes through the international date line * may return a bad bounding box.//from w w w . ja va 2s. co m * * @param polyAry array of coordinates (lon,lat,lon,lat,lon,lat..etc) * @return Array of bounding box coordinates in the following order: West South East North. Also * described as minX, minY, maxX, maxY (where longitude is the X-axis, and latitude is the * Y-axis). */ private double[] createBBoxFromPolygon(String[] polyAry) { double minX = Double.POSITIVE_INFINITY; double minY = Double.POSITIVE_INFINITY; double maxX = Double.NEGATIVE_INFINITY; double maxY = Double.NEGATIVE_INFINITY; double curX, curY; for (int i = 0; i < polyAry.length - 1; i += 2) { LOGGER.debug("polyToBBox: lon - {} lat - {}", polyAry[i], polyAry[i + 1]); curX = Double.parseDouble(polyAry[i]); curY = Double.parseDouble(polyAry[i + 1]); if (curX < minX) { minX = curX; } if (curX > maxX) { maxX = curX; } if (curY < minY) { minY = curY; } if (curY > maxY) { maxY = curY; } } return new double[] { minX, minY, maxX, maxY }; }
From source file:dkpro.similarity.algorithms.wikipedia.measures.WikipediaSimilarityMeasureBase.java
/** * Relatedness values < 0 are invalid by definition. * @param valueList A list containing the values. * @return Returns the minimum of the values in the given list, or NO_SENSE if the value list is empty. *///ww w. java 2 s . com public double getMinimum(List<Double> valueList) { if (valueList == null) { return NO_SENSE; } boolean valid = false; double minimum = Double.POSITIVE_INFINITY; for (double value : valueList) { if (value >= 0 && value < minimum) { minimum = value; valid = true; } } if (!valid) { return NO_SENSE; } return minimum; }
From source file:com.genericworkflownodes.knime.config.writer.CTDConfigurationWriter.java
private void addDoubleListParameterRestrictions(Parameter<?> p, StringBuffer restriction) { DoubleListParameter dlp = (DoubleListParameter) p; boolean lbSet = Double.NEGATIVE_INFINITY != dlp.getLowerBound().doubleValue(); boolean ubSet = Double.POSITIVE_INFINITY != dlp.getUpperBound().doubleValue(); if (lbSet) {/* ww w .j ava 2 s . c o m*/ restriction.append(String.format(Locale.ENGLISH, "%f", dlp.getLowerBound()) .replaceAll(REMOVE_TRAILING_0_RE, "").replaceAll(REMOVE_TRAILING_DOT_RE, "")); } if (ubSet || lbSet) { restriction.append(':'); } if (ubSet) { restriction.append(String.format(Locale.ENGLISH, "%f", dlp.getUpperBound()) .replaceAll(REMOVE_TRAILING_0_RE, "").replaceAll(REMOVE_TRAILING_DOT_RE, "")); } }
From source file:com.opengamma.analytics.financial.model.finitedifference.CraigSneydFiniteDifference2D.java
private int sor(final int steps, final double[][] v, final double[] q, final double[][] mx, final int j, final double omega) { double sum;/*from w w w . ja va 2s. co m*/ int min; int max; int count = 0; double scale = 1.0; double errorSqr = Double.POSITIVE_INFINITY; while (errorSqr / (scale + 1e-10) > 1e-18 && count < 1000) { errorSqr = 0.0; scale = 0.0; for (int l = 0; l <= steps; l++) { min = (l == steps ? 0 : Math.max(0, l - 1)); max = (l == 0 ? steps : Math.min(steps, 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] * v[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; v[l][j] += correction; scale += v[l][j] * v[l][j]; } count++; } return count; }
From source file:org.jax.maanova.fit.gui.ResidualPlotPanel.java
private void mouseMoved(MouseEvent e) { if (this.showTooltip) { Point2D chartPoint = this.chartPanel.toChartPoint(e.getPoint()); // find the nearest probe XYProbeData[] xyProbeData = this.getXYData(); double nearestDistance = Double.POSITIVE_INFINITY; int nearestArrayIndex = -1; int nearestDotIndex = -1; for (int arrayIndex = 0; arrayIndex < xyProbeData.length; arrayIndex++) { double[] currXData = xyProbeData[arrayIndex].getXData(); double[] currYData = xyProbeData[arrayIndex].getYData(); for (int dotIndex = 0; dotIndex < currXData.length; dotIndex++) { double currDist = chartPoint.distanceSq(currXData[dotIndex], currYData[dotIndex]); if (currDist < nearestDistance) { nearestDistance = currDist; nearestArrayIndex = arrayIndex; nearestDotIndex = dotIndex; }/*ww w . j av a 2 s .c o m*/ } } if (nearestArrayIndex == -1) { this.clearProbePopup(); } else { XYProbeData nearestArrayData = xyProbeData[nearestArrayIndex]; Point2D probeJava2DCoord = this.getJava2DCoordinates(nearestArrayData.getXData()[nearestDotIndex], nearestArrayData.getYData()[nearestDotIndex]); double java2DDist = probeJava2DCoord.distance(e.getX(), e.getY()); // is the probe close enough to be worth showing (in pixel distance) if (java2DDist <= PlotUtil.SCATTER_PLOT_DOT_SIZE_PIXELS * 2) { this.showProbePopup(nearestArrayIndex, nearestArrayData.getProbeIndices()[nearestDotIndex], nearestArrayData.getXData()[nearestDotIndex], nearestArrayData.getYData()[nearestDotIndex], e.getX(), e.getY()); } else { this.clearProbePopup(); } } } }
From source file:com.addthis.hydra.data.tree.prop.DataReservoir.java
private ValueObject generateValueObject(String key, String separator) { long targetEpoch = -1; int numObservations = -1; double sigma = Double.POSITIVE_INFINITY; double percentile = 0; boolean doubleToLongBits = false; int minMeasurement = Integer.MIN_VALUE; boolean raw = false; String mode = "sigma"; if (key == null) { return null; }/*w w w .j ava 2 s .c o m*/ String[] kvpairs = key.split(Pattern.quote("~")); for (String kvpair : kvpairs) { String[] kv = kvpair.split(Pattern.quote(separator)); if (kv.length == 2) { String kvkey = kv[0]; String kvvalue = kv[1]; switch (kvkey) { case "double": doubleToLongBits = Boolean.parseBoolean(kvvalue); break; case "epoch": targetEpoch = Long.parseLong(kvvalue); break; case "sigma": sigma = Double.parseDouble(kvvalue); break; case "min": minMeasurement = Integer.parseInt(kvvalue); break; case "obs": numObservations = Integer.parseInt(kvvalue); break; case "raw": raw = Boolean.parseBoolean(kvvalue); break; case "percentile": percentile = Double.parseDouble(kvvalue); break; case "mode": mode = kvvalue; break; default: throw new RuntimeException("Unknown key " + kvkey); } } } if (mode.equals("get")) { long count = retrieveCount(targetEpoch); if (count < 0) { return ValueFactory.create(0L); } else { return ValueFactory.create(count); } } else { DataReservoirValue.Builder builder = new DataReservoirValue.Builder(); builder.setTargetEpoch(targetEpoch); builder.setNumObservations(numObservations); builder.setDoubleToLongBits(doubleToLongBits); builder.setRaw(raw); builder.setSigma(sigma); builder.setMinMeasurement(minMeasurement); builder.setPercentile(percentile); builder.setMode(mode); DataReservoirValue value = builder.build(this); return value; } }
From source file:com.vsthost.rnd.jdeoptim.DE.java
/** * Returns the index of the minimum element in the double array provided. * * @param array The array to be search for the index of the minimum element. * @return The index of the minimum element in the array. *///from w w w. j a va 2 s . c o m private static int minIndex(double[] array) { // Declare and initialize the return value: int index = -1; // Declare and initialize the minimum value: double min = Double.POSITIVE_INFINITY; // Find the best score and update values: for (int i = 0; i < array.length; i++) { if (array[i] <= min) { min = array[i]; index = i; } } // Done, return: return index; }
From source file:clus.statistic.ClassificationStat.java
/** * Currently only used to compute the default dispersion within rule heuristics. *//* w w w. ja va 2 s . c om*/ public double getDispersion(ClusAttributeWeights scale, RowData data) { System.err.println(getClass().getName() + ": getDispersion(): Not yet implemented!"); return Double.POSITIVE_INFINITY; }
From source file:edu.duke.cs.osprey.ematrix.epic.SeriesFitter.java
static double[] fitSeriesIterative(DoubleMatrix1D[] samp, double trueVals[], double weights[], double lambda, boolean includeConst, int order, double bCutoffs[], double bCutoffs2[], int PCOrder, boolean isPC[]) { long startTime = System.currentTimeMillis(); System.out.println("Starting fitSeriesIterative..."); int numSamples = samp.length; int nd = samp[0].size(); if (bCutoffs.length == 1) {//single bCutoff to be used double bCutoff = bCutoffs[0]; bCutoffs = new double[numSamples]; Arrays.fill(bCutoffs, bCutoff); }//w w w . ja va 2s. c om //now bCutoffs has a cutoff for each sample int numParams = getNumParams(nd, includeConst, order); if (PCOrder > order) {//add in parameters for PC orders int numPCs = countTrue(isPC); for (int n = order + 1; n <= PCOrder; n++) numParams += getNumParamsForOrder(numPCs, n); } //now set up the data for the iterative fits //samples can be turned on and off using fitWeights //we will need to create two entries for samples with trueVals between bCutoff and bCutoff2 //since they may be turned on either to penalize deviation from bCutoff or from trueVal //first entry for each entry will penalize deviation from bCutoff if trueVal>=bCutoff //secondEntries are for trueVals between bCutoff and bCutoff2 ArrayList<Integer> secondEntries = new ArrayList<>();//list of samples needing second entry HashMap<Integer, Integer> revSecondEntries = new HashMap<>();//reverse lookup for (int s = 0; s < numSamples; s++) { if ((trueVals[s] >= bCutoffs[s]) && (trueVals[s] < bCutoffs2[s])) {//as in isRestraintActive revSecondEntries.put(s, secondEntries.size()); secondEntries.add(s); } } int numRestraints = numSamples + secondEntries.size(); //data for basic least-squares fits DoubleMatrix1D[] fitSamp = new DoubleMatrix1D[numRestraints]; double fitTrueVals[] = new double[numRestraints]; double fitWeights[] = new double[numRestraints]; for (int s = 0; s < numSamples; s++) {//"normal" entries fitSamp[s] = samp[s]; if (trueVals[s] >= bCutoffs[s]) { fitWeights[s] = 0; fitTrueVals[s] = bCutoffs[s]; } else { fitWeights[s] = weights[s]; fitTrueVals[s] = trueVals[s]; } } for (int s2 = 0; s2 < secondEntries.size(); s2++) { fitSamp[numSamples + s2] = samp[secondEntries.get(s2)]; fitWeights[numSamples + s2] = 0; fitTrueVals[numSamples + s2] = trueVals[secondEntries.get(s2)]; } //Initial guess of set P is all points with trueVals[s] >= bCutoff //that is, all points that have possible series values that make the restraint inactive boolean done = false; double coeffs[] = null; double meanResidual = 0, weightSum = 0; double prevResid = Double.POSITIVE_INFINITY; double oldCoeffs[] = null; double oldSerVals[] = new double[numSamples];//values of series at each sample, for previous iteration //preallocating to all infinity because all trueVals[s]>=bCutoff points start outside P Arrays.fill(oldSerVals, Double.POSITIVE_INFINITY); //for updating boolean firstFit = true;//first fit is not an update DoubleMatrix1D c = DoubleFactory1D.dense.make(numParams);//matrices we update (used in fit) DoubleMatrix2D M = DoubleFactory2D.dense.make(numParams, numParams); double oldFitWeights[] = null; //double fitWeightsCheck[] = fitWeights.clone();//DEBUG!!! while (!done) { if (firstFit) { coeffs = fitSeries(fitSamp, fitTrueVals, fitWeights, lambda, includeConst, order, PCOrder, isPC, false, c, M); firstFit = false; } else { double weightDiffs[] = fitWeights.clone(); for (int s = 0; s < numRestraints; s++) weightDiffs[s] -= oldFitWeights[s]; coeffs = fitSeries(fitSamp, fitTrueVals, weightDiffs, lambda, includeConst, order, PCOrder, isPC, true, c, M); //DEBUG!!! /* for(int s=0; s<numRestraints; s++){ fitWeightsCheck[s] += weightDiffs[s]; if(fitWeightsCheck[s] != fitWeights[s]){ int cefAO = 1111; } } double checkCoeffs[] = fitSeries(fitSamp, fitTrueVals, fitWeights, lambda, includeConst, order, PCOrder, isPC, false, null, null); for(int a=0; a<coeffs.length; a++){ if(Math.abs(checkCoeffs[a]-coeffs[a])>1e-10){ int abc=123; } } DoubleMatrix1D c2 = DoubleFactory1D.dense.make(numParams); DoubleMatrix2D M2 = DoubleFactory2D.dense.make(numParams,numParams); double checkCoeffs2[] = fitSeries(fitSamp, fitTrueVals, fitWeights, lambda, includeConst, order, PCOrder, isPC, false, c2, M2); for(int a=0; a<coeffs.length; a++){ if(Math.abs(checkCoeffs2[a]-coeffs[a])>1e-10){ int abc=123; } }*/ //DEBUG!!! } oldFitWeights = fitWeights.clone(); done = true; ArrayList<SampleCutoffCrossing> scc = new ArrayList<SampleCutoffCrossing>(); meanResidual = 0; weightSum = 0; //boolean doneNoTol = true, done2 = true;//DEBUG!!! //values of series at each sample, based on coeffs double serVals[] = new double[numSamples]; for (int s = 0; s < numSamples; s++) { serVals[s] = evalSeries(coeffs, samp[s], nd, includeConst, order, PCOrder, isPC); //If each series value is below or above bcutoff according to whether the //coeffs were generated by fitting with or without that value included //(respectively), then we have found a local and thus the global minimum //in this quadratic piece of the objective function, so we're done //check for doneness first, using tolerance //i.e. are coeffs (derived from fitWeights) consistent with fitWeights, //within numerical error? If so we have a global minimum if (trueVals[s] >= bCutoffs[s]) { if (fitWeights[s] > 0) { if (!isRestraintTypeActive(trueVals[s], serVals[s], bCutoffs[s], bCutoffs2[s], false, -1e-6)) done = false;//fitWeights penalizing deviation from bCutoff, and this isn't right at coeffs } else { boolean secondRestraintOn = revSecondEntries.containsKey(s); if (secondRestraintOn) secondRestraintOn = (fitWeights[numSamples + revSecondEntries.get(s)] > 0); if (secondRestraintOn) { if (!isRestraintTypeActive(trueVals[s], serVals[s], bCutoffs[s], bCutoffs2[s], true, -1e-6)) done = false;//fitWeights penalizing deviation from trueVal, and this isn't right at coeffs } else {//restraints currently off if (isRestraintActive(trueVals[s], serVals[s], bCutoffs[s], bCutoffs2[s], 1e-6)) done = false;//a restraint should be on at coeffs } } } //for trueVals below bCutoff, restraints don't turn on and off //DEBUG!!!!! //trying to calculate done w/o tol /* if( trueVals[s]>=bCutoffs[s] ){ if(fitWeights[s]>0){ if(!isRestraintTypeActive(trueVals[s],serVals[s],bCutoffs[s],bCutoffs2[s],false,0)) doneNoTol = false;//fitWeights penalizing deviation from bCutoff, and this isn't right at coeffs } else { boolean secondRestraintOn = revSecondEntries.containsKey(s); if(secondRestraintOn) secondRestraintOn = (fitWeights[numSamples+revSecondEntries.get(s)]>0); if(secondRestraintOn){ if(!isRestraintTypeActive(trueVals[s],serVals[s],bCutoffs[s],bCutoffs2[s],true,0)) doneNoTol = false;//fitWeights penalizing deviation from trueVal, and this isn't right at coeffs } else {//restraints currently off if(isRestraintActive(trueVals[s],serVals[s],bCutoffs[s],bCutoffs2[s],0)) doneNoTol = false;//a restraint should be on at coeffs } } } //OK now done2 will be calculated and should be the same but will be calculated like //the weight changes below if( trueVals[s]>=bCutoffs[s] ){ if(isRestraintTypeActive(trueVals[s],serVals[s],bCutoffs[s],bCutoffs2[s],false)){ //activate penalty for deviating from bCutoff if(fitWeights[s]!=weights[s]) done2 = false; if(revSecondEntries.containsKey(s)){ if(fitWeights[numSamples+revSecondEntries.get(s)]!=0) done2 = false; } } else if(isRestraintTypeActive(trueVals[s],serVals[s],bCutoffs[s],bCutoffs2[s],true)){ //activate penalty for deviating from trueVal if(fitWeights[s] != 0) done2 = false; if(revSecondEntries.containsKey(s)){ if(fitWeights[numSamples+revSecondEntries.get(s)] != weights[s]) done2 = false; } else throw new RuntimeException("ERROR: should have second entry for restraint but don't!!"); } else { //deactivate all penalties if(fitWeights[s] != 0) done2 = false; if(revSecondEntries.containsKey(s)){ if(fitWeights[numSamples+revSecondEntries.get(s)] != 0) done2 = false; } //no contribution to residual } }*/ //DEBUG!!! //Now calculate mean residual and crossing points, and update fitWeights double residTerm = 0; if (trueVals[s] >= bCutoffs[s]) { if (isRestraintTypeActive(trueVals[s], serVals[s], bCutoffs[s], bCutoffs2[s], false)) { //activate penalty for deviating from bCutoff fitWeights[s] = weights[s]; if (revSecondEntries.containsKey(s)) fitWeights[numSamples + revSecondEntries.get(s)] = 0; residTerm = (serVals[s] - bCutoffs[s]) * (serVals[s] - bCutoffs[s]); } else if (isRestraintTypeActive(trueVals[s], serVals[s], bCutoffs[s], bCutoffs2[s], true)) { //activate penalty for deviating from trueVal fitWeights[s] = 0; if (revSecondEntries.containsKey(s)) fitWeights[numSamples + revSecondEntries.get(s)] = weights[s]; else throw new RuntimeException("ERROR: should have second entry for restraint but don't!!"); residTerm = (serVals[s] - trueVals[s]) * (serVals[s] - trueVals[s]); } else { //deactivate all penalties fitWeights[s] = 0; if (revSecondEntries.containsKey(s)) fitWeights[numSamples + revSecondEntries.get(s)] = 0; //no contribution to residual } } else //normal least-squares penalty. fitWeights[s] will stay at weights[s] residTerm = (serVals[s] - trueVals[s]) * (serVals[s] - trueVals[s]); meanResidual += weights[s] * residTerm; //If want sample-by-sample output... //System.out.println("TRAININGSET TRUE: "+trueVals[s]+" SER: "+serVals[s]); weightSum += weights[s]; } meanResidual /= weightSum; if (meanResidual == prevResid) System.out.println(); //DEBUG!!! /* if(done!=doneNoTol || done2!=done){ //Let's see what happens if we remove the tolerance... done = doneNoTol; } if(done){ double checkCoeffs[] = fitSeries(fitSamp, fitTrueVals, fitWeights, lambda, includeConst, order, PCOrder, isPC, false, null, null); for(int a=0; a<coeffs.length; a++){ if(Math.abs(checkCoeffs[a]-coeffs[a])>1e-10){ int abc=123; } } }*/ //DEBUG!!! if ((!done) && (meanResidual >= prevResid)) { //Did not obtain a decrease using the Newton step //Let's do an exact line search to rectify the situation if (!useLineSearch) { System.out.println("Skipping line search, returning with residual " + prevResid); return oldCoeffs; } System.out.println("LINE SEARCH"); for (int s = 0; s < numSamples; s++) { //If we go in or out of either type of restraint between serVals and oldSerVals, we create //a SampleCutoffCrossing of the appropriate type (upper or lower (ordinary) restraint) if ((isRestraintTypeActive(trueVals[s], oldSerVals[s], bCutoffs[s], bCutoffs2[s], false)) != (isRestraintTypeActive(trueVals[s], serVals[s], bCutoffs[s], bCutoffs2[s], false))) { //If the restraint disappears at one end we know trueVal>=bCutoff here //create lower restraint SampleCutoffCrossing double crossingPoint = (bCutoffs[s] - oldSerVals[s]) / (serVals[s] - oldSerVals[s]); scc.add(new SampleCutoffCrossing(s, crossingPoint, false)); } if ((isRestraintTypeActive(trueVals[s], oldSerVals[s], bCutoffs[s], bCutoffs2[s], true)) != (isRestraintTypeActive(trueVals[s], serVals[s], bCutoffs[s], bCutoffs2[s], true))) { //upper double crossingPoint2 = (trueVals[s] - oldSerVals[s]) / (serVals[s] - oldSerVals[s]); scc.add(new SampleCutoffCrossing(s, crossingPoint2, true)); } } int changeCount = scc.size(); Collections.sort(scc); TreeSet<Integer> crossingIndices = new TreeSet<Integer>(); for (SampleCutoffCrossing cr : scc) { int s = cr.sampleIndex; crossingIndices.add(s); if (cr.upperRestraint) {//penalizing difference from trueVal cr.quadTerm = weights[s] * (serVals[s] - oldSerVals[s]) * (serVals[s] - oldSerVals[s]); cr.linTerm = 2 * weights[s] * (serVals[s] - oldSerVals[s]) * (oldSerVals[s] - trueVals[s]); cr.constTerm = weights[s] * (oldSerVals[s] - trueVals[s]) * (oldSerVals[s] - trueVals[s]); } else {//penalizing difference from lesser of bCutoff, trueVal cr.quadTerm = weights[s] * (serVals[s] - oldSerVals[s]) * (serVals[s] - oldSerVals[s]); double baseVal = Math.min(trueVals[s], bCutoffs[s]); cr.linTerm = 2 * weights[s] * (serVals[s] - oldSerVals[s]) * (oldSerVals[s] - baseVal); cr.constTerm = weights[s] * (oldSerVals[s] - baseVal) * (oldSerVals[s] - baseVal); //cr.linTerm = 2*weights[s]*(serVals[s]-oldSerVals[s])*(oldSerVals[s]-trueVals[s]); //cr.constTerm = weights[s]*(oldSerVals[s]-trueVals[s])*(oldSerVals[s]-trueVals[s]); } } //Set up quadratic function double quadTerm = 0; double linTerm = 0; double constTerm = 0; //Add in contributions from all non-cutoff-crossing points for (int s = 0; s < numSamples; s++) { if (!crossingIndices.contains(s)) { if (isRestraintTypeActive(trueVals[s], oldSerVals[s], bCutoffs[s], bCutoffs2[s], false)) { //if(trueVals[s]<bCutoffs[s]||serVals[s]<bCutoffs[s]){//penalizing difference from lesser of bCutoff, trueVal quadTerm += weights[s] * (serVals[s] - oldSerVals[s]) * (serVals[s] - oldSerVals[s]); double baseVal = Math.min(trueVals[s], bCutoffs[s]); linTerm += 2 * weights[s] * (serVals[s] - oldSerVals[s]) * (oldSerVals[s] - baseVal); constTerm += weights[s] * (oldSerVals[s] - baseVal) * (oldSerVals[s] - baseVal); //linTerm += 2*weights[s]*(serVals[s]-oldSerVals[s])*(oldSerVals[s]-trueVals[s]); //constTerm += weights[s]*(oldSerVals[s]-trueVals[s])*(oldSerVals[s]-trueVals[s]); } else if (isRestraintTypeActive(trueVals[s], oldSerVals[s], bCutoffs[s], bCutoffs2[s], true)) { //else if(isRestraintActive(trueVals[s],serVals[s],bCutoffs[s],bCutoffs2[s],0)){//penalizing difference from trueVal quadTerm += weights[s] * (serVals[s] - oldSerVals[s]) * (serVals[s] - oldSerVals[s]); linTerm += 2 * weights[s] * (serVals[s] - oldSerVals[s]) * (oldSerVals[s] - trueVals[s]); constTerm += weights[s] * (oldSerVals[s] - trueVals[s]) * (oldSerVals[s] - trueVals[s]); } } } //contributions from cutoff-crossing points at the beginning of the interval //(i.e. at coeffs) for (SampleCutoffCrossing cr : scc) { int s = cr.sampleIndex; if (isRestraintTypeActive(trueVals[s], serVals[s], bCutoffs[s], bCutoffs2[s], cr.upperRestraint)) { //check if this particular restraint (rather than either restraint for this s) is active //if(serVals[s]<bCutoffs[s]){ quadTerm += cr.quadTerm; linTerm += cr.linTerm; constTerm += cr.constTerm; } } //double checkMeanResid = (quadTerm+linTerm+constTerm)/weightSum;//evaluate objective function at a=1 //should match meanResidual! double prevNodeResid = Double.POSITIVE_INFINITY; //The first increase we may consider is from node 0 to 1 //Now walk back until we get an increase //then the minimum will be in one of the last two quadratic pieces int lowestNodeIndex = 0; for (int curChangeIndex = changeCount - 1; curChangeIndex >= 0; curChangeIndex--) { SampleCutoffCrossing cr = scc.get(curChangeIndex); double a = cr.crossingPoint; double curNodeResid = (a * a * quadTerm + a * linTerm + constTerm) / weightSum; int s = cr.sampleIndex; if (curNodeResid > prevNodeResid) { lowestNodeIndex = curChangeIndex + 1; break; } //if the restraint is being removed going from serVals to oldSerVals, remove its coefficients from the quadratic function //if it's being added, add them if (isRestraintTypeActive(trueVals[s], serVals[s], bCutoffs[s], bCutoffs2[s], cr.upperRestraint)) { quadTerm -= cr.quadTerm; linTerm -= cr.linTerm; constTerm -= cr.constTerm; } else { quadTerm += cr.quadTerm; linTerm += cr.linTerm; constTerm += cr.constTerm; } prevNodeResid = curNodeResid; } //At this point, we know our minimum is in either the piece with the //current quad, lin, constTerms or the one we looked at right before //(where lowestNodeIndex is the node separating them) double a_min = -linTerm / (2 * quadTerm); SampleCutoffCrossing cr = scc.get(lowestNodeIndex); if (a_min > cr.crossingPoint) {//minimum must be in previous piece //revert quad and linTerms int s = cr.sampleIndex; //change back quadratic-function coefficients if (isRestraintTypeActive(trueVals[s], serVals[s], bCutoffs[s], bCutoffs2[s], cr.upperRestraint)) { quadTerm += cr.quadTerm; linTerm += cr.linTerm; constTerm += cr.constTerm; } else { quadTerm -= cr.quadTerm; linTerm -= cr.linTerm; constTerm -= cr.constTerm; } a_min = -linTerm / (2 * quadTerm); } //double minResid = (a_min*a_min*quadTerm + a_min*linTerm + constTerm)/weightSum; for (int p = 0; p < numParams; p++) coeffs[p] = coeffs[p] * a_min + (1 - a_min) * oldCoeffs[p]; double minResid = 0; for (int s = 0; s < numSamples; s++) { serVals[s] = evalSeries(coeffs, samp[s], nd, includeConst, order, PCOrder, isPC); if (trueVals[s] >= bCutoffs[s]) { if (isRestraintTypeActive(trueVals[s], serVals[s], bCutoffs[s], bCutoffs2[s], false)) { //activate penalty for deviating from bCutoff fitWeights[s] = weights[s]; if (revSecondEntries.containsKey(s)) fitWeights[numSamples + revSecondEntries.get(s)] = 0; minResid += weights[s] * (serVals[s] - bCutoffs[s]) * (serVals[s] - bCutoffs[s]); } else if (isRestraintTypeActive(trueVals[s], serVals[s], bCutoffs[s], bCutoffs2[s], true)) { minResid += weights[s] * (serVals[s] - trueVals[s]) * (serVals[s] - trueVals[s]); //activate penalty for deviating from trueVal fitWeights[s] = 0; fitWeights[numSamples + revSecondEntries.get(s)] = weights[s]; } else { //deactivate all penalties fitWeights[s] = 0; if (revSecondEntries.containsKey(s)) fitWeights[numSamples + revSecondEntries.get(s)] = 0; } } else minResid += weights[s] * (serVals[s] - trueVals[s]) * (serVals[s] - trueVals[s]); } minResid /= weightSum; if (minResid >= prevResid) { //consider to have converged //NOTE THIS CAN HAPPEN IF THE QUADRATIC APPROXIMATION AT OLDCOEFFS HAS SOLUTION //FAR FROM THE EXACT VALUE (ASSUMING EXACT FITSERIES) OF 1 //THIS CAN HAPPEN IF WE'RE GETTING BELOW THE NUMERICAL PRECISION OF FITSERIES System.out.println("TRAINING SET MEAN RESIDUAL:" + prevResid); System.out.println("CONVERGED IN LINE SEARCH, line search min: " + minResid); System.out.println("fitSeriesIterative time (ms): " + (System.currentTimeMillis() - startTime)); return oldCoeffs; } meanResidual = minResid; } /* //trying backtracking line search, hoping it'll be more numerically stable double backtrackRate = 0.5; double a = 1;//how far we want to be on the line from oldCoeffs (which should allow some descent //along the line) to coeffs (which overshoots) double minResid; double aCoeffs[] = new double[coeffs.length];//coefficients at our point on the line do { a *= backtrackRate; for(int p=0; p<numParams; p++) aCoeffs[p] = coeffs[p]*a + (1-a)*oldCoeffs[p]; minResid = 0; for(int s=0; s<numSamples; s++){ serVals[s] = evalSeries(aCoeffs,samp[s],nd,includeConst,order,PCOrder,isPC); if( trueVals[s]>=bCutoffs[s] ){ if(isRestraintTypeActive(trueVals[s],serVals[s],bCutoffs[s],bCutoffs2[s],false)){ minResid += weights[s]*(serVals[s]-bCutoffs[s])*(serVals[s]-bCutoffs[s]); weights2[s] = weights[s]; } else if(isRestraintTypeActive(trueVals[s],serVals[s],bCutoffs[s],bCutoffs2[s],true)){ minResid += weights[s]*(serVals[s]-trueVals[s])*(serVals[s]-trueVals[s]); weights2[s] = weights[s]; } else weights2[s] = 0; } else minResid += weights[s]*(serVals[s]-trueVals[s])*(serVals[s]-trueVals[s]); } minResid /= weightSum; } while(minResid>prevResid); if(a<1e-4) System.out.println("Warning: line search a got down to "+a); meanResidual = minResid; } */ oldCoeffs = coeffs; System.out.println("STEP RESIDUAL: " + meanResidual); prevResid = meanResidual; oldSerVals = serVals; } System.out.println("TRAINING SET MEAN RESIDUAL:" + meanResidual); System.out.println("fitSeriesIterative time (ms): " + (System.currentTimeMillis() - startTime)); //DEBUG!!! //Note this assumes weights are all 1! /* int numCoeffs = coeffs.length; DoubleMatrix1D residGrad = DoubleFactory1D.dense.make(numCoeffs); double checkResid = 0; for(int s=0; s<samp.length; s++){ double diff = 0; double fitVal = evalSeries(coeffs,samp[s],nd,includeConst,order,PCOrder,isPC); if(trueVals[s]>=bCutoffs[s] && fitVal<bCutoffs[s]){ diff = fitVal-bCutoffs[s]; if(fitWeights[s]==0){ int allosaurus = 15; } if(revSecondEntries.containsKey(s)){ if(fitWeights[numSamples+revSecondEntries.get(s)]>0){ int allosaurus=15; } } } else if( trueVals[s]<bCutoffs[s] || (trueVals[s]<bCutoffs2[s] && fitVal>trueVals[s]) ){ diff = fitVal-trueVals[s]; if(trueVals[s]<bCutoffs[s]){ if(fitWeights[s]==0){ int allosaurus = 15; } } else if(fitWeights[numSamples+revSecondEntries.get(s)]==0 || fitWeights[s]>0){ int allosaurus=15; } } else if(fitWeights[s]>0){ int allosaurus=15; } else if(revSecondEntries.containsKey(s)){ if(fitWeights[numSamples+revSecondEntries.get(s)]>0){ int allosaurus=15; } } checkResid += diff*diff; DoubleMatrix1D termMonomials = DoubleFactory1D.dense.make(numCoeffs); SeriesFitter.calcSampParamCoeffs( termMonomials, samp[s], nd, includeConst, order, PCOrder, isPC ); residGrad.assign( termMonomials, Functions.plusMult(2*diff) ); } checkResid /= samp.length; residGrad.assign( Functions.mult(1./samp.length) ); double checkCoeffs[] = fitSeries(fitSamp, fitTrueVals, fitWeights, lambda, includeConst, order, PCOrder, isPC, false, null, null); double checkResid2 = 0; DoubleMatrix1D residGrad2 = DoubleFactory1D.dense.make(numCoeffs); for(int s=0; s<samp.length; s++){ double diff = 0; double fitVal = evalSeries(checkCoeffs,samp[s],nd,includeConst,order,PCOrder,isPC); if(trueVals[s]>=bCutoffs[s] && fitVal<bCutoffs[s]){ diff = fitVal-bCutoffs[s]; if(fitWeights[s]==0){ int allosaurus = 15; } if(revSecondEntries.containsKey(s)){ if(fitWeights[numSamples+revSecondEntries.get(s)]>0){ int allosaurus=15; } } } else if( trueVals[s]<bCutoffs[s] || (trueVals[s]<bCutoffs2[s] && fitVal>trueVals[s]) ){ diff = fitVal-trueVals[s]; if(trueVals[s]<bCutoffs[s]){ if(fitWeights[s]==0){ int allosaurus = 15; } } else if(fitWeights[numSamples+revSecondEntries.get(s)]==0 || fitWeights[s]>0){ int allosaurus=15; } } else if(fitWeights[s]>0){ int allosaurus=15; } else if(revSecondEntries.containsKey(s)){ if(fitWeights[numSamples+revSecondEntries.get(s)]>0){ int allosaurus=15; } } checkResid2 += diff*diff; DoubleMatrix1D termMonomials = DoubleFactory1D.dense.make(numCoeffs); SeriesFitter.calcSampParamCoeffs( termMonomials, samp[s], nd, includeConst, order, PCOrder, isPC ); residGrad2.assign( termMonomials, Functions.plusMult(2*diff) ); } checkResid2 /= samp.length; residGrad2.assign( Functions.mult(1./samp.length) ); if(residGrad.zDotProduct(residGrad)>1e-10){ int struthiomimus = 23; } int cubist = -1; */ //DEBUG!! return coeffs; }