List of usage examples for java.util Arrays binarySearch
public static int binarySearch(Object[] a, int fromIndex, int toIndex, Object key)
From source file:io.warp10.continuum.gts.GTSOutliersHelper.java
/** * Applying Seasonal Entropy Hybrid test * This test is based on piecewise decomposition where trend components are approximated by median and seasonal components by entropy of the cycle sub-series. * An ESD test is passed upon the residuals. * /*from w w w . j a v a2s . c o m*/ * It differs from hybridTest by approximating seasonal component instead of using STL. * But in many cases this approximation is more useful than estimation of STL. * * @param gts * @param k Upper bound of suspected number of outliers * @param alpha Significance level with which to accept or reject anomalies. Default is 0.05 * * @return anomalous_ticks * * @throws WarpScriptException */ public static List<Long> entropyHybridTest(GeoTimeSerie gts, int buckets_per_period, int periods_per_piece, int k, double alpha) throws WarpScriptException { doubleCheck(gts); List<Long> anomalous_ticks = new ArrayList<Long>(); if (!GTSHelper.isBucketized(gts)) { throw new WarpScriptException("GTS must be bucketized"); } if (k >= periods_per_piece * buckets_per_period / 2) { throw new WarpScriptException( "Upper bound of number of outliers must be less than half of the number of observations per piece"); } // // SubSerie attributes // GeoTimeSerie subgts = null; GeoTimeSerie subsubgts = null; GeoTimeSerie seasonal = null; // number of pieces long pieces = gts.bucketcount / buckets_per_period / periods_per_piece; // number of buckets per piece int bpp = periods_per_piece * buckets_per_period; long lb = gts.lastbucket; long bs = gts.bucketspan; for (int u = 0; u < pieces; u++) { long start = lb - bs * ((pieces - u) * bpp - 1); long stop = lb - bs * (pieces - u - 1) * bpp; // we don't start from the first bucket subgts = GTSHelper.subSerie(gts, start, stop, false, false, subgts); subgts.lastbucket = stop; subgts.bucketcount = bpp; subgts.bucketspan = bs; // entropy seasonal extraction if (null == seasonal) { seasonal = new GeoTimeSerie(bpp); seasonal.doubleValues = new double[bpp]; seasonal.ticks = new long[bpp]; } else { GTSHelper.reset(seasonal); } seasonal.type = TYPE.DOUBLE; // TODO(JCV): make this a method ? for (int v = 0; v < buckets_per_period; v++) { subsubgts = GTSHelper.subCycleSerie(subgts, stop - v * bs, buckets_per_period, true, subsubgts); // compute entropy of absolute modified zscore double[] madsigma = madsigma(subsubgts, true); double median = madsigma[0]; double mad = madsigma[1]; double sum = 0.0D; for (int w = 0; w < subsubgts.values; w++) { subsubgts.doubleValues[w] = 0.0D != mad ? Math.abs((subsubgts.doubleValues[w] - median) / mad) : 1.0D; sum += subsubgts.doubleValues[w]; } double entropy = 0.0D; for (int w = 0; w < subsubgts.values; w++) { subsubgts.doubleValues[w] /= sum; double tmp = subsubgts.doubleValues[w]; if (0.0D != tmp) { entropy -= tmp * Math.log(tmp); } } // normalize entropy and handle case where all items are the same if (0.0D != entropy) { entropy /= Math.log(subsubgts.values); } else { entropy = 1.0D; } // update seasonal for (int w = 0; w < subsubgts.values; w++) { GTSHelper.setValue(seasonal, subsubgts.ticks[w], entropy * subsubgts.doubleValues[w]); } } GTSHelper.sort(seasonal); double m = median(seasonal); int idx = 0; for (int i = 0; i < subgts.values; i++) { idx = Arrays.binarySearch(seasonal.ticks, idx, seasonal.values, subgts.ticks[i]); if (idx < 0) { throw new WarpScriptException("Internal bug method entropyHybridTest: can't find tick " + subgts.ticks[i] + " in seasonal.ticks"); } else { subgts.doubleValues[i] -= (seasonal.doubleValues[idx] + m); } } anomalous_ticks.addAll(ESDTest(subgts, k, true, alpha)); } return anomalous_ticks; }
From source file:org.csploit.android.services.UpdateService.java
/** * extract an archive into a directory//from w w w. jav a 2 s. c o m * * @throws IOException if some I/O error occurs * @throws java.util.concurrent.CancellationException if task is cancelled by user * @throws java.lang.InterruptedException when the the running thread get cancelled. */ private void extract() throws RuntimeException, IOException, InterruptedException, ChildManager.ChildNotStartedException { ArchiveInputStream is = null; ArchiveEntry entry; CountingInputStream counter; OutputStream outputStream = null; File f, inFile; File[] list; String name; String envPath; final StringBuffer sb = new StringBuffer(); int mode; int count; long total; boolean isTar, r, w, x, isElf, isScript; short percentage, old_percentage; Child which; if (mCurrentTask.path == null || mCurrentTask.outputDir == null) return; mBuilder.setContentTitle(getString(R.string.extracting)).setContentText("").setContentInfo("") .setSmallIcon(android.R.drawable.ic_popup_sync).setProgress(100, 0, false); mNotificationManager.notify(NOTIFICATION_ID, mBuilder.build()); Logger.info(String.format("extracting '%s' to '%s'", mCurrentTask.path, mCurrentTask.outputDir)); envPath = null; which = null; try { if (mCurrentTask.fixShebang) { which = System.getTools().raw.async("which env", new Raw.RawReceiver() { @Override public void onNewLine(String line) { sb.delete(0, sb.length()); sb.append(line); } }); } inFile = new File(mCurrentTask.path); total = inFile.length(); counter = new CountingInputStream(new FileInputStream(inFile)); is = openArchiveStream(counter); isTar = mCurrentTask.archiver.equals(archiveAlgorithm.tar); old_percentage = -1; f = new File(mCurrentTask.outputDir); if (f.exists() && f.isDirectory() && (list = f.listFiles()) != null && list.length > 2) wipe(); if (mCurrentTask.fixShebang) { if (execShell(which, "cancelled while retrieving env path") != 0) { throw new RuntimeException("cannot find 'env' executable"); } envPath = sb.toString(); } while (mRunning && (entry = is.getNextEntry()) != null) { name = entry.getName().replaceFirst("^\\./?", ""); if (mCurrentTask.skipRoot) { if (name.contains("/")) name = name.substring(name.indexOf('/') + 1); else if (entry.isDirectory()) continue; } f = new File(mCurrentTask.outputDir, name); isElf = isScript = false; if (entry.isDirectory()) { if (!f.exists()) { if (!f.mkdirs()) { throw new IOException( String.format("Couldn't create directory '%s'.", f.getAbsolutePath())); } } } else { byte[] buffer = null; byte[] writeMe = null; outputStream = new FileOutputStream(f); // check il file is an ELF or a script if ((!isTar || mCurrentTask.fixShebang) && entry.getSize() > 4) { writeMe = buffer = new byte[4]; IOUtils.readFully(is, buffer); if (buffer[0] == 0x7F && buffer[1] == 0x45 && buffer[2] == 0x4C && buffer[3] == 0x46) { isElf = true; } else if (buffer[0] == '#' && buffer[1] == '!') { isScript = true; ByteArrayOutputStream firstLine = new ByteArrayOutputStream(); int newline = -1; // assume that '\n' is more far then 4 chars. firstLine.write(buffer); buffer = new byte[1024]; count = 0; while (mRunning && (count = is.read(buffer)) >= 0 && (newline = Arrays.binarySearch(buffer, 0, count, (byte) 0x0A)) < 0) { firstLine.write(buffer, 0, count); } if (!mRunning) { throw new CancellationException("cancelled while searching for newline."); } else if (count < 0) { newline = count = 0; } else if (newline < 0) { newline = count; } firstLine.write(buffer, 0, newline); firstLine.close(); byte[] newFirstLine = new String(firstLine.toByteArray()) .replace("/usr/bin/env", envPath).getBytes(); writeMe = new byte[newFirstLine.length + (count - newline)]; java.lang.System.arraycopy(newFirstLine, 0, writeMe, 0, newFirstLine.length); java.lang.System.arraycopy(buffer, newline, writeMe, newFirstLine.length, count - newline); } } if (writeMe != null) { outputStream.write(writeMe); } IOUtils.copy(is, outputStream); outputStream.close(); outputStream = null; percentage = (short) (((double) counter.getBytesRead() / total) * 100); if (percentage != old_percentage) { mBuilder.setProgress(100, percentage, false).setContentInfo(percentage + "%"); mNotificationManager.notify(NOTIFICATION_ID, mBuilder.build()); old_percentage = percentage; } } // Zip does not store file permissions. if (isTar) { mode = ((TarArchiveEntry) entry).getMode(); r = (mode & 0400) > 0; w = (mode & 0200) > 0; x = (mode & 0100) > 0; } else if (isElf || isScript) { r = w = x = true; } else { continue; } if (!f.setExecutable(x, true)) { Logger.warning(String.format("cannot set executable permission of '%s'", name)); } if (!f.setWritable(w, true)) { Logger.warning(String.format("cannot set writable permission of '%s'", name)); } if (!f.setReadable(r, true)) { Logger.warning(String.format("cannot set readable permission of '%s'", name)); } } if (!mRunning) throw new CancellationException("extraction cancelled."); Logger.info("extraction completed"); f = new File(mCurrentTask.outputDir, ".nomedia"); if (f.createNewFile()) Logger.info(".nomedia created"); mBuilder.setContentInfo("").setProgress(100, 100, true); mNotificationManager.notify(NOTIFICATION_ID, mBuilder.build()); } finally { if (is != null) is.close(); if (outputStream != null) outputStream.close(); } }
From source file:io.warp10.continuum.gts.GTSHelper.java
public static int indexAtTick(GeoTimeSerie gts, long tick) { if (0 == gts.values) { return -1; }/*from ww w. j a va 2 s . c o m*/ sort(gts, false); // // Attempt to locate the tick // int idx = Arrays.binarySearch(gts.ticks, 0, gts.values, tick); if (idx < 0) { return -1; } return idx; }
From source file:io.warp10.continuum.gts.GTSHelper.java
/** * Return the value in a Geo Time Serie at a given timestamp. * //from ww w .jav a2 s . c o m * The GeoTimeSerie instance will be sorted if it is not already * * @param gts GeoTimeSerie instance from which to extract value * @param tick Timestamp at which to read the value * @return The value at 'tick' or null if none exists */ public static Object valueAtTick(GeoTimeSerie gts, long tick) { if (0 == gts.values) { return null; } // // Force sorting in natural ordering of ticks // sort(gts, false); // // Attempt to locate the tick // int idx = Arrays.binarySearch(gts.ticks, 0, gts.values, tick); if (idx < 0) { return null; } else { if (TYPE.LONG == gts.type) { return gts.longValues[idx]; } else if (TYPE.DOUBLE == gts.type) { return gts.doubleValues[idx]; } else if (TYPE.STRING == gts.type) { return gts.stringValues[idx]; } else if (TYPE.BOOLEAN == gts.type) { return gts.booleanValues.get(idx); } else { return null; } } }
From source file:io.warp10.continuum.gts.GTSHelper.java
/** * Return the location in a Geo Time Serie at a given timestamp. * // w w w . j av a 2s.co m * The GeoTimeSerie instance will be sorted if it is not already. * * @param gts GeoTimeSerie instance from which to extract location * @param tick Timestamp at which to read the location * @return The location at 'tick' (NO_LOCATION if none set). */ public static long locationAtTick(GeoTimeSerie gts, long tick) { if (null == gts.locations) { return GeoTimeSerie.NO_LOCATION; } sort(gts, false); int idx = Arrays.binarySearch(gts.ticks, 0, gts.values, tick); if (idx < 0) { return GeoTimeSerie.NO_LOCATION; } else { return gts.locations[idx]; } }
From source file:io.warp10.continuum.gts.GTSHelper.java
/** * Return the elevation in a Geo Time Serie at a given timestamp. * //from www .j a va 2 s .c o m * The GeoTimeSerie instance will be sorted if it is not already. * * @param gts GeoTimeSerie instance from which to extract elevation * @param tick Timestamp at which to read the elevation * @return The elevation at 'tick' (NO_ELEVATION if none set). */ public static long elevationAtTick(GeoTimeSerie gts, long tick) { if (null == gts.elevations) { return GeoTimeSerie.NO_ELEVATION; } sort(gts, false); int idx = Arrays.binarySearch(gts.ticks, 0, gts.values, tick); if (idx < 0) { return GeoTimeSerie.NO_ELEVATION; } else { return gts.elevations[idx]; } }
From source file:org.csploit.android.core.UpdateService.java
/** * extract an archive into a directory//from w ww .jav a 2s. com * * @throws IOException if some I/O error occurs * @throws java.util.concurrent.CancellationException if task is cancelled by user * @throws java.lang.InterruptedException when the the running thread get cancelled. */ private void extract() throws CancellationException, RuntimeException, IOException, InterruptedException, ChildManager.ChildNotStartedException { ArchiveInputStream is = null; ArchiveEntry entry; CountingInputStream counter; OutputStream outputStream = null; File f, inFile; File[] list; String name; String envPath; final StringBuffer sb = new StringBuffer(); int mode; int count; long total; boolean isTar, r, w, x, isElf, isScript; short percentage, old_percentage; Child which; DiffMatchPatch dmp; if (mCurrentTask.path == null || mCurrentTask.outputDir == null) return; mBuilder.setContentTitle(getString(R.string.extracting)).setContentText("").setContentInfo("") .setSmallIcon(android.R.drawable.ic_popup_sync).setProgress(100, 0, false); mNotificationManager.notify(NOTIFICATION_ID, mBuilder.build()); Logger.info(String.format("extracting '%s' to '%s'", mCurrentTask.path, mCurrentTask.outputDir)); envPath = null; which = null; try { if (mCurrentTask.fixShebang) { which = System.getTools().raw.async("which env", new Raw.RawReceiver() { @Override public void onNewLine(String line) { sb.delete(0, sb.length()); sb.append(line); } }); } inFile = new File(mCurrentTask.path); total = inFile.length(); counter = new CountingInputStream(new FileInputStream(inFile)); is = openArchiveStream(counter); isTar = mCurrentTask.archiver.equals(archiveAlgorithm.tar); old_percentage = -1; dmp = (mCurrentTask.patches != null && mCurrentTask.patches.size() > 0) ? new DiffMatchPatch() : null; f = new File(mCurrentTask.outputDir); if (f.exists() && f.isDirectory() && (list = f.listFiles()) != null && list.length > 2) wipe(); if (mCurrentTask.fixShebang) { if (execShell(which, "cancelled while retrieving env path") != 0) { throw new RuntimeException("cannot find 'env' executable"); } envPath = sb.toString(); } while (mRunning && (entry = is.getNextEntry()) != null) { name = entry.getName().replaceFirst("^\\./?", ""); if (mCurrentTask.skipRoot) { if (name.contains("/")) name = name.substring(name.indexOf('/') + 1); else if (entry.isDirectory()) continue; } f = new File(mCurrentTask.outputDir, name); isElf = isScript = false; if (entry.isDirectory()) { if (!f.exists()) { if (!f.mkdirs()) { throw new IOException( String.format("Couldn't create directory '%s'.", f.getAbsolutePath())); } } } else { byte[] buffer = null; byte[] writeMe = null; // patch the file if (dmp != null && mCurrentTask.patches.containsKey(name)) { buffer = new byte[(int) entry.getSize()]; IOUtils.readFully(is, buffer); writeMe = buffer = ((String) dmp.patch_apply(mCurrentTask.patches.get(name), new String(buffer))[0]).getBytes(); } outputStream = new FileOutputStream(f); // check il file is an ELF or a script if ((!isTar || mCurrentTask.fixShebang) && entry.getSize() > 4) { if (buffer == null) { writeMe = buffer = new byte[4]; IOUtils.readFully(is, buffer); if (buffer[0] == 0x7F && buffer[1] == 0x45 && buffer[2] == 0x4C && buffer[3] == 0x46) { isElf = true; } else if (buffer[0] == '#' && buffer[1] == '!') { isScript = true; ByteArrayOutputStream firstLine = new ByteArrayOutputStream(); int newline = -1; // assume that '\n' is more far then 4 chars. firstLine.write(buffer); buffer = new byte[1024]; count = 0; while (mRunning && (count = is.read(buffer)) >= 0 && (newline = Arrays.binarySearch(buffer, 0, count, (byte) 0x0A)) < 0) { firstLine.write(buffer, 0, count); } if (!mRunning) { throw new CancellationException("cancelled while searching for newline."); } else if (count < 0) { newline = count = 0; } else if (newline < 0) { newline = count; } firstLine.write(buffer, 0, newline); firstLine.close(); byte[] newFirstLine = new String(firstLine.toByteArray()) .replace("/usr/bin/env", envPath).getBytes(); writeMe = new byte[newFirstLine.length + (count - newline)]; java.lang.System.arraycopy(newFirstLine, 0, writeMe, 0, newFirstLine.length); java.lang.System.arraycopy(buffer, newline, writeMe, newFirstLine.length, count - newline); } } else { if (buffer[0] == 0x7F && buffer[1] == 0x45 && buffer[2] == 0x4C && buffer[3] == 0x46) { isElf = true; } else if (buffer[0] == '#' && buffer[1] == '!') { isScript = true; int newline = Arrays.binarySearch(buffer, (byte) 0x0A); if (newline < 0) newline = buffer.length; byte[] newFirstLine = new String(buffer, 0, newline) .replace("/usr/bin/env", envPath).getBytes(); writeMe = new byte[buffer.length + (newFirstLine.length - newline)]; java.lang.System.arraycopy(newFirstLine, 0, writeMe, 0, newFirstLine.length); java.lang.System.arraycopy(buffer, newline, writeMe, newFirstLine.length, newFirstLine.length - newline); } } } if (writeMe != null) { outputStream.write(writeMe); } IOUtils.copy(is, outputStream); outputStream.close(); outputStream = null; percentage = (short) (((double) counter.getBytesRead() / total) * 100); if (percentage != old_percentage) { mBuilder.setProgress(100, percentage, false).setContentInfo(percentage + "%"); mNotificationManager.notify(NOTIFICATION_ID, mBuilder.build()); old_percentage = percentage; } } // Zip does not store file permissions. if (isTar) { mode = ((TarArchiveEntry) entry).getMode(); r = (mode & 0400) > 0; w = (mode & 0200) > 0; x = (mode & 0100) > 0; } else if (isElf || isScript) { r = w = x = true; } else { continue; } if (!f.setExecutable(x, true)) { Logger.warning(String.format("cannot set executable permission of '%s'", name)); } if (!f.setWritable(w, true)) { Logger.warning(String.format("cannot set writable permission of '%s'", name)); } if (!f.setReadable(r, true)) { Logger.warning(String.format("cannot set readable permission of '%s'", name)); } } if (!mRunning) throw new CancellationException("extraction cancelled."); Logger.info("extraction completed"); f = new File(mCurrentTask.outputDir, ".nomedia"); if (f.createNewFile()) Logger.info(".nomedia created"); mBuilder.setContentInfo("").setProgress(100, 100, true); mNotificationManager.notify(NOTIFICATION_ID, mBuilder.build()); } finally { if (is != null) is.close(); if (outputStream != null) outputStream.close(); } }
From source file:io.warp10.continuum.gts.GTSHelper.java
/** * Return a new GeoTimeSerie instance containing only the value of 'gts' * which fall between 'timestamp' (inclusive) and 'timestamp' - 'span' (exclusive) * //w w w .j a v a 2 s .c o m * The resulting GTS instance will be sorted. * * @param gts GeoTimeSerie from which to extract values. * @param starttimestamp Oldest timestamp to consider (in microseconds) * @param stoptimestamp Most recent timestamp to consider (in microseconds) * @param overwrite Should we overwrite measurements which occur at the same timestamp to only keep the last one added. * @param copyLabels If true, labels will be copied from the original GTS to the subserie. * * @return The computed sub Geo Time Serie */ public static final GeoTimeSerie subSerie(GeoTimeSerie gts, long starttimestamp, long stoptimestamp, boolean overwrite, boolean copyLabels, GeoTimeSerie subgts) { // FIXME(hbs): should a subserie of a bucketized GTS be also bucketized? With possible non existant values. This would impact Mappers/Bucketizers // // Create sub GTS // // The size hint impacts performance, choose it wisely... if (null == subgts) { subgts = new GeoTimeSerie(128); // // Copy name and labels // subgts.setName(gts.getName()); if (copyLabels) { subgts.setLabels(gts.getLabels()); } } else { GTSHelper.reset(subgts); } if (null == gts.ticks || 0 == gts.values) { return subgts; } // // Sort GTS so ticks are ordered // GTSHelper.sort(gts); // // Determine index to start at // int lastidx = Arrays.binarySearch(gts.ticks, 0, gts.values, stoptimestamp); // // The upper timestamp is less than the first tick, so subserie is necessarly empty // if (-1 == lastidx) { return subgts; } else if (lastidx < 0) { lastidx = -lastidx - 1; if (lastidx >= gts.values) { lastidx = gts.values - 1; } } int firstidx = Arrays.binarySearch(gts.ticks, 0, lastidx + 1, starttimestamp); if (firstidx < 0) { firstidx = -firstidx - 1; } if (firstidx >= gts.values) { return subgts; } // // Extract values/locations/elevations that lie in the requested interval // for (int i = firstidx; i <= lastidx; i++) { if (gts.ticks[i] >= starttimestamp && gts.ticks[i] <= stoptimestamp) { setValue(subgts, gts.ticks[i], null != gts.locations ? gts.locations[i] : GeoTimeSerie.NO_LOCATION, null != gts.elevations ? gts.elevations[i] : GeoTimeSerie.NO_ELEVATION, valueAtIndex(gts, i), overwrite); } } return subgts; }
From source file:io.warp10.continuum.gts.GTSHelper.java
/** * Return a new GeoTimeSerie instance containing only the values of 'gts' * that fall under timestamps that are bound by the same modulo class. * //from w w w . jav a 2s. c o m * The resulting GTS instance will be sorted and bucketized. * * @param gts GeoTimeSerie from which to extract values. It must be bucketized. * @param lastbucket Most recent timestamp to consider (in microseconds) * @param buckets_per_period Number of buckets of the input gts that sum up to a bucket for the sub cycle gts * @param overwrite Should we overwrite measurements which occur at the same timestamp to only keep the last one added. * * @return The computed sub cycle Geo Time Serie */ public static final GeoTimeSerie subCycleSerie(GeoTimeSerie gts, long lastbucket, int buckets_per_period, boolean overwrite, GeoTimeSerie subgts) throws WarpScriptException { if (!isBucketized(gts)) { throw new WarpScriptException("GTS must be bucketized"); } if (0 != (gts.lastbucket - lastbucket) % gts.bucketspan) { throw new WarpScriptException( "lasbucket parameter of subCycleSerie method must fall on an actual bucket of the gts input"); } // // Create sub GTS // // The size hint impacts performance, choose it wisely... if (null == subgts) { subgts = new GeoTimeSerie(lastbucket, (gts.bucketcount - (int) ((gts.lastbucket - lastbucket) / gts.bucketspan) - 1) / buckets_per_period + 1, gts.bucketspan * buckets_per_period, (int) Math.max(1.4 * gts.bucketcount, gts.sizehint) / buckets_per_period); } else { subgts.values = 0; subgts.type = TYPE.UNDEFINED; subgts.lastbucket = lastbucket; subgts.bucketcount = (gts.bucketcount - (int) ((gts.lastbucket - lastbucket) / gts.bucketspan) - 1) / buckets_per_period + 1; subgts.bucketspan = gts.bucketspan * buckets_per_period; } if (null == gts.ticks || 0 == gts.values) { return subgts; } // // For each tick, search if tick is in gts, then copy it to subgts, else noop // Iterator<Long> iter = tickIterator(subgts, true); long tick; sort(gts); int i = gts.values; while (iter.hasNext()) { tick = iter.next(); i = Arrays.binarySearch(gts.ticks, 0, i, tick); if (i >= 0) { setValue(subgts, gts.ticks[i], null != gts.locations ? gts.locations[i] : GeoTimeSerie.NO_LOCATION, null != gts.elevations ? gts.elevations[i] : GeoTimeSerie.NO_ELEVATION, valueAtIndex(gts, i), overwrite); } } return subgts; }
From source file:io.warp10.continuum.gts.GTSHelper.java
/** * Compute fast and robust version of LOWESS on a Geo Time Series, * with a polynomial fit of degree p > 0. * /*from w w w . j a v a 2 s. c om*/ * @see http://www.stat.washington.edu/courses/stat527/s14/readings/Cleveland_JASA_1979.pdf * * @param gts Input GTS * @param q Bandwith, i.e. number of nearest neighbours to consider when applying LOWESS * @param r Robustness, i.e. number of robustifying iterations * @param d Delta in s, i.e. acceptable neighbourhood radius within which LOWESS is not recomputed * close points are approximated by polynomial interpolation, d should remain < 0.1*(lasttick-firstick) in most cases * @param p Degree, i.e. degree of the polynomial fit * best usage p=1 or p=2 ; it will likely return an overflow error if p is too big * * @param weights : optional array that store the weights * @param rho : optional array that store the robustness weights * @param inplace : should the gts returned be the same object than the input * * @return a smoothed GTS * @throws WarpScriptException */ public static GeoTimeSerie rlowess(GeoTimeSerie gts, int q, int r, long d, int p, double[] weights, double[] rho, boolean inplace) throws WarpScriptException { if (TYPE.DOUBLE != gts.type && TYPE.LONG != gts.type) { throw new WarpScriptException("Can only smooth numeric Geo Time Series."); } if (q < 1) { throw new WarpScriptException("Bandwitdth parameter must be greater than 0"); } if (r < 0) { throw new WarpScriptException("Robustness parameter must be greater or equal to 0"); } if (d < 0) { throw new WarpScriptException("Delta parameter must be greater or equal to 0"); } if (p < 1) { throw new WarpScriptException("Degree of polynomial fit must be greater than 0"); } if (p > 9) { throw new WarpScriptException("Degree of polynomial fit should remain small (lower than 10)"); } // // Sort the ticks // sort(gts, false); // // Check if number of missing values is reasonable // Note that (bucketcount - values) is not the number of missing values but a minorant // if (gts.bucketcount - gts.values > 500000) { throw new WarpScriptException("More than 500000 missing values"); } // // Check for duplicate ticks // long previous = -1; for (int t = 0; t < gts.values; t++) { long current = gts.ticks[t]; if (previous == current) { throw new WarpScriptException("Can't be applied on GTS with duplicate ticks"); } previous = current; } // // At each step of robustifying operations, // we compute values in transient_smoothed using updated weights rho * weights // int size = isBucketized(gts) ? gts.bucketcount : gts.values; int sizehint = Math.max(gts.sizehint, (int) 1.1 * size); double[] transient_smoothed = new double[sizehint]; int nvalues = q < size ? q : size; // Allocate an array for weights if (null == weights) { weights = new double[nvalues]; } else { if (weights.length < nvalues) { throw new WarpScriptException("in rlowess weights array too small"); } } // Allocate an array to update the weights through the robustness iterations and another for the absolute of the residuals // The weights used will be rho*weights double[] residual; if (r > 0) { if (null == rho) { rho = new double[gts.values]; Arrays.fill(rho, 1.0D); } else { if (rho.length < nvalues) { throw new WarpScriptException("in rlowess rho array too small"); } } residual = new double[gts.values]; } else { residual = null; } // Regression parameters double[] beta = new double[p + 1]; // // Robustifying iterations // int r_iter = 0; while (r_iter < r + 1) { // // In order to speed up the computations, // we will skip some points pointed by iter. // We use iter_follower to interpolate the skipped points afterward. // Iterator<Long> iter = tickIterator(gts, false); Iterator<Long> iter_follower; if (0.0 == d) { iter_follower = null; } else { iter_follower = tickIterator(gts, false); } // Index in the ticks/values array of the input gts int idx = 0; // Index in the ticks/values array of the output gts (values are also estimated for null points of the input) int ridx = 0; int ridx_last = 0; // Last tick estimated (set to -d-1 so (tick-last)>d at first iter) and its index in the result long last = d * (-1) - 1; int idx_last = 0; // // When we find a tick that is not within distance d of the last estimated tick, // then either we estimate it, // or if at least one tick has been skipped just before, // then we estimate the last skipped one and interpolate the others. // We then take back the loop from the former last skipped tick. // long last_skipped = 0; boolean skip = false; // Have skipped ticks been interpolated in last loop ? boolean resolved = false; // Current tick long tick = 0; while (iter.hasNext() || resolved) { if (!resolved) { tick = iter.next(); } else { resolved = false; } // Skip points that are too close from the previous estimated one, unless its the last if (iter.hasNext() && (tick - last <= d)) { last_skipped = tick; skip = true; ridx++; } else { if (!skip) { // advance idx to the first neighbour at the right whose value is not null while (idx < gts.values - 1 && tick > tickAtIndex(gts, idx)) { idx++; } // compute value at tick transient_smoothed[ridx] = pointwise_lowess(gts, idx, tick, nvalues, p, weights, rho, beta); // update residual if tick had a non-null value if (r_iter < r && tick == tickAtIndex(gts, idx)) { residual[idx] = Math.abs( ((Number) valueAtIndex(gts, idx)).doubleValue() - transient_smoothed[ridx]); } if (null != iter_follower) { iter_follower.next(); last = tick; idx_last = idx; ridx_last = ridx; } ridx++; } else { if (!iter.hasNext() && (tick - last <= d)) { last_skipped = tick; ridx++; } // advance idx to the first neighbour at the right whose value is not null while (idx < gts.values - 1 && last_skipped > tickAtIndex(gts, idx)) { idx++; } // compute value at last_skipped tick transient_smoothed[ridx - 1] = pointwise_lowess(gts, idx, last_skipped, nvalues, p, weights, rho, beta); // update residual if tick had a non-null value if (r_iter < r && last_skipped == tickAtIndex(gts, idx)) { residual[idx] = Math.abs( ((Number) valueAtIndex(gts, idx)).doubleValue() - transient_smoothed[ridx - 1]); } // // Linear interpolation of skipped points // double denom = last_skipped - last; long skipped = iter_follower.next(); int ridx_s = ridx_last + 1; while (last_skipped > skipped) { // interpolate double alpha = (skipped - last) / denom; transient_smoothed[ridx_s] = alpha * transient_smoothed[ridx - 1] + (1 - alpha) * transient_smoothed[ridx_last]; // update residual if tick had a non-null value int sidx; if (r_iter < r && 0 < (sidx = Arrays.binarySearch(gts.ticks, idx_last, idx, skipped))) { residual[sidx] = Math.abs(((Number) valueAtIndex(gts, sidx)).doubleValue() - transient_smoothed[ridx_s]); } skipped = iter_follower.next(); ridx_s++; } if (iter.hasNext() || (tick - last > d)) { //updates skip = false; resolved = true; last = last_skipped; idx_last = idx; ridx_last = ridx - 1; } } } } // // Update robustifying weights (except last time or if r is 0) // if (r_iter < r) { // rho's values will be recomputed anyway so sorted can borrow its body double[] sorted = rho; sorted = Arrays.copyOf(residual, gts.values); Arrays.sort(sorted); // compute median of abs(residual) double median; if (gts.values % 2 == 0) { median = (sorted[gts.values / 2] + sorted[gts.values / 2 - 1]) / 2; } else { median = sorted[gts.values / 2]; } // compute h = 6 * median and rho = bisquare(|residual|/h) double h = 6 * median; for (int k = 0; k < gts.values; k++) { if (0 == h) { rho[k] = 1.0D; } else { double u = residual[k] / h; if (u >= 1.0) { rho[k] = 0.0D; } else { rho[k] = 1.0D - u * u; rho[k] = rho[k] * rho[k]; } } } } r_iter++; } // // Copying result to output // boolean hasLocations = null != gts.locations; boolean hasElevations = null != gts.elevations; // // We separate case without or with missing values. // if (!isBucketized(gts) || (gts.values == gts.bucketcount && gts.lastbucket == gts.ticks[gts.values - 1] && gts.lastbucket - gts.bucketspan * (gts.bucketcount - 1) == gts.ticks[0])) { if (inplace) { if (TYPE.LONG == gts.type) { gts.longValues = null; gts.type = TYPE.DOUBLE; } gts.doubleValues = transient_smoothed; return gts; } else { GeoTimeSerie smoothed = gts.cloneEmpty(sizehint); try { smoothed.reset(Arrays.copyOf(gts.ticks, sizehint), hasLocations ? Arrays.copyOf(gts.locations, sizehint) : null, hasElevations ? Arrays.copyOf(gts.elevations, sizehint) : null, transient_smoothed, size); } catch (IOException ioe) { throw new WarpScriptException("IOException in reset method: " + ioe.getMessage()); } return smoothed; } } else { // // Case with missing values // if (inplace) { // TODO(JCV): need to unit test location and elevation settings in this case if (hasLocations && gts.locations.length != sizehint) { gts.locations = Arrays.copyOf(gts.locations, sizehint); } if (hasElevations && gts.elevations.length != sizehint) { gts.elevations = Arrays.copyOf(gts.elevations, sizehint); } if (gts.ticks.length != sizehint) { gts.ticks = Arrays.copyOf(gts.ticks, sizehint); } // We try to allocate the lesser additionnal memory so we // fill locations and elevations backward with values upfront in the same array if (hasLocations || hasElevations) { Iterator<Long> iter = tickIterator(gts, true); int idx = gts.values - 1; int idx_new = gts.bucketcount; while (iter.hasNext()) { long tick = iter.next(); idx_new--; // search if tick was present in original gts while (idx > 0 && tick < gts.ticks[idx]) { idx--; } if (hasLocations) { gts.locations[idx_new] = tick == gts.ticks[idx] ? gts.locations[idx] : GeoTimeSerie.NO_LOCATION; } if (hasElevations) { gts.elevations[idx_new] = tick == gts.ticks[idx] ? gts.locations[idx] : GeoTimeSerie.NO_ELEVATION; } gts.ticks[idx_new] = tick; } } if (TYPE.LONG == gts.type) { gts.longValues = null; gts.type = TYPE.DOUBLE; } gts.doubleValues = transient_smoothed; gts.values = size; return gts; } else { GeoTimeSerie smoothed = gts.cloneEmpty(sizehint); // ticks long[] ticks = new long[sizehint]; int idx = 0; Iterator<Long> iter = tickIterator(gts, false); while (iter.hasNext()) { ticks[idx] = iter.next(); idx++; } // locations int v = 0; long[] locations = null; if (hasLocations) { locations = new long[sizehint]; for (int u = 0; u < size; u++) { v = Arrays.binarySearch(gts.ticks, v, gts.values, smoothed.ticks[u]); locations[u] = v < 0 ? GeoTimeSerie.NO_LOCATION : gts.locations[v]; } } // elevations v = 0; long[] elevations = null; if (hasElevations) { elevations = new long[sizehint]; for (int u = 0; u < size; u++) { v = Arrays.binarySearch(gts.ticks, v, gts.values, smoothed.ticks[u]); elevations[u] = v < 0 ? GeoTimeSerie.NO_ELEVATION : gts.elevations[v]; } } try { smoothed.reset(ticks, locations, elevations, transient_smoothed, size); } catch (IOException ioe) { throw new WarpScriptException("IOException in reset method: " + ioe.getMessage()); } return smoothed; } } }