Skip to content

Add stable history detection using ROC test #39

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
nhey wants to merge 2 commits into
base: roc
Choose a base branch
from
Open

Add stable history detection using ROC test #39

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
3f44d21
Convert to f64 and add stable history detection using ROC test.
May 7, 2021
34a2c11
Generated Futhark OpenCL code.
May 7, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions bfast/base.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,7 @@ class BFASTMonitorBase(ABC):
def __init__(
self,
start_monitor,
history="all",
freq=365,
k=3,
hfrac=0.25,
Expand All @@ -21,6 +22,7 @@ def __init__(
verbose=0,
):
self.start_monitor = start_monitor
self.history = history
self.freq = freq
self.k = k
self.hfrac = hfrac
Expand Down
34 changes: 33 additions & 1 deletion bfast/models.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,14 @@ class BFASTMonitor():
A datetime object specifying the start of
the monitoring phase.

history : str, default "all"
Specification of the start of a stable history period.
Can be one of "ROC" and "all". If set to "ROC", a
reverse-ordered CUSUM test will be employed to
automatically detect the start of a stable history
period. If set to "all", the start of the stable
history period is the beginning of the time series.

freq : int, default 365
The frequency for the seasonal model.

Expand Down Expand Up @@ -89,6 +97,7 @@ class BFASTMonitor():
def __init__(
self,
start_monitor,
history="all",
freq=365,
k=3,
hfrac=0.25,
Expand All @@ -102,7 +111,11 @@ def __init__(
detailed_results=False,
find_magnitudes=True,
):
history_enum = ("all", "ROC")
if history not in history_enum:
raise Exception("Current implementation can only handle the following values for history: {}".format(history_enum))
self.start_monitor = start_monitor
self.history = history
self.freq = freq
self.k = k
self.hfrac = hfrac
Expand Down Expand Up @@ -143,6 +156,7 @@ def fit(self, data, dates, n_chunks=None, nan_value=0):
if self.backend == 'python':
self._model = BFASTMonitorPython(
start_monitor=self.start_monitor,
history=self.history,
freq=self.freq,
k=self.k,
hfrac=self.hfrac,
Expand All @@ -155,6 +169,7 @@ def fit(self, data, dates, n_chunks=None, nan_value=0):
elif self.backend == 'python-mp':
self._model = BFASTMonitorPython(
start_monitor=self.start_monitor,
history=self.history,
freq=self.freq,
k=self.k,
hfrac=self.hfrac,
Expand All @@ -168,6 +183,7 @@ def fit(self, data, dates, n_chunks=None, nan_value=0):
elif self.backend == 'opencl':
self._model = BFASTMonitorOpenCL(
start_monitor=self.start_monitor,
history=self.history,
freq=self.freq,
k=self.k,
hfrac=self.hfrac,
Expand Down Expand Up @@ -312,9 +328,25 @@ def valids(self):
-------
array-like : An array containing the number
of valid values for each pixel in the
aray data
array data
"""
if self._is_fitted():
return self._model.valids

raise Exception("Model not yet fitted!")

@property
def history_starts(self):
""" Returns index of the start of the stable
history period for each pixel

Returns
-------
array-like : An array containing indices
of starts of stable history periods
in the array data
"""
if self._is_fitted():
return self._model.history_starts

raise Exception("Model not yet fitted!")
75 changes: 55 additions & 20 deletions bfast/monitor/opencl/base.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,9 @@
import pyopencl.array as pycl_array

from bfast.base import BFASTMonitorBase
from bfast.monitor.utils import compute_end_history, compute_lam, map_indices
from bfast.monitor.utils import (
compute_end_history, compute_lam, map_indices, compute_lam_brownian
)

# generated by Futhark
from .bfastfinal import bfastfinal
Expand All @@ -37,6 +39,13 @@ class BFASTMonitorOpenCL(BFASTMonitorBase):
start_monitor : datetime
A datetime object specifying the start of
the monitoring phase.
history : str, default "all"
Specification of the start of a stable history period.
Can be one of "ROC" and "all". If set to "ROC", a
reverse-ordered CUSUM test will be employed to
automatically detect the start of a stable history
period. If set to "all", the start of the stable
history period is the beginning of the time series.
freq : int, default 365
The frequency for the seasonal model
k : int, default 3
Expand Down Expand Up @@ -71,6 +80,7 @@ class BFASTMonitorOpenCL(BFASTMonitorBase):
"""
def __init__(self,
start_monitor,
history="all",
freq=365,
k=3,
hfrac=0.25,
Expand All @@ -89,6 +99,7 @@ def __init__(self,
raise Exception("Current implementation can only handle the following values for k: {}".format(k_valid))

super().__init__(start_monitor,
history,
freq,
k=k,
hfrac=hfrac,
Expand All @@ -102,14 +113,17 @@ def __init__(self,
self.platform_id = platform_id
self.device_id = device_id

history_enum = { "ROC": -1, "all": 0 }
self.history_num = history_enum[self.history]

# initialize device
self._init_device(platform_id, device_id)

self.futobj = bfastfinal(command_queue=self.queue,
interactive=False,
default_tile_size=8,
default_reg_tile_size=3,
sizes=self._get_futhark_params())
default_reg_tile_size=3)
# sizes=self._get_futhark_params())

def _init_device(self, platform_id, device_id):
""" Initializes the device."""
Expand Down Expand Up @@ -224,6 +238,7 @@ def fit(self, data, dates, n_chunks=None, nan_value=0):
self.breaks = results['breaks']
self.means = results['means']
self.valids = results['valids']
self.history_starts = results['history_starts']

if self.find_magnitudes or self.detailed_results:
self.magnitudes = results['magnitudes']
Expand Down Expand Up @@ -291,10 +306,11 @@ def _fit_single(self, data, dates, nan_value):

def _fit_single_preprocess(self, data, dates, nan_value):
start = time.time()
mapped_indices = map_indices(dates).astype(numpy.int32)
mapped_indices = map_indices(dates).astype(numpy.int64)
N = data.shape[0]
self.n = compute_end_history(dates, self.start_monitor)
self.lam = compute_lam(N, self.hfrac, self.level, self.period)
self.conf_ROC = compute_lam_brownian(self.level)
end = time.time()

if self.verbose > 0:
Expand Down Expand Up @@ -339,32 +355,43 @@ def _fit_single_kernel(self, y_cl, mapped_indices_cl):
means, \
magnitudes, \
y_error, \
y_pred = self.futobj.mainDetailed(trend,
self.k,
self.n,
self.freq,
self.hfrac,
self.lam,
mapped_indices_cl, y_cl)
y_pred, \
hist = self.futobj.mainDetailed(trend,
self.k,
self.n,
self.freq,
self.hfrac,
self.level,
self.lam,
self.history_num,
self.conf_ROC,
mapped_indices_cl, y_cl)
elif self.find_magnitudes:
Ns, \
breaks, \
means, \
magnitudes = self.futobj.mainMagnitude(trend,
self.k,
self.n,
self.freq,
self.hfrac,
self.lam,
mapped_indices_cl, y_cl)
else:
Ns, breaks, means = self.futobj.main(trend,
magnitudes, \
hist = self.futobj.mainMagnitude(trend,
self.k,
self.n,
self.freq,
self.hfrac,
self.level,
self.lam,
self.history_num,
self.conf_ROC,
mapped_indices_cl, y_cl)
else:
Ns, breaks, means, hist = self.futobj.main(trend,
self.k,
self.n,
self.freq,
self.hfrac,
self.level,
self.lam,
self.history_num,
self.conf_ROC,
mapped_indices_cl, y_cl)

end = time.time()

Expand All @@ -380,6 +407,7 @@ def _fit_single_kernel(self, y_cl, mapped_indices_cl):
results['breaks'] = breaks
results['means'] = means
results['valids'] = Ns
results['history_starts'] = hist

if self.find_magnitudes or self.detailed_results:
results['magnitudes'] = magnitudes
Expand All @@ -398,6 +426,7 @@ def _fit_single_postprocessing(self, results, oshape):
results['breaks'] = results['breaks'].get().reshape(oshape[1:])
results['means'] = results['means'].get().reshape(oshape[1:])
results['valids'] = results['valids'].get().T.reshape(oshape[1:])
results['history_starts'] = results['history_starts'].get().T.reshape(oshape[1:])

if self.find_magnitudes or self.detailed_results:
results['magnitudes'] = results['magnitudes'].get().reshape(oshape[1:])
Expand Down Expand Up @@ -435,6 +464,12 @@ def __append_results(self, res, results):
else:
results['valids'] = res['valids']

if 'history_starts' in results.keys():
results['history_starts'] = numpy.concatenate([results['history_starts'],
res['history_starts']], axis=0)
else:
results['history_starts'] = res['history_starts']

if self.find_magnitudes or self.detailed_results:
if 'magnitudes' in results.keys():
results['magnitudes'] = numpy.concatenate([results['magnitudes'], res['magnitudes']], axis=0)
Expand Down
Loading