Skip to content

Commit

Permalink
Convert to f64 and add stable history detection using ROC test.
Browse files Browse the repository at this point in the history
  • Loading branch information
Nikolaj Hey Hinnerskov authored and Nikolaj Hey Hinnerskov committed Jul 13, 2021
1 parent d4c1511 commit a37c899
Show file tree
Hide file tree
Showing 15 changed files with 761 additions and 196 deletions.
2 changes: 2 additions & 0 deletions bfast/base.py
View file
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
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
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
4 changes: 2 additions & 2 deletions bfast/monitor/opencl/futhark/Makefile
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
bfastfinal.py: bfastfinal.fut helpers.fut
futhark pkg add github.com/diku-dk/sorts
bfastfinal.py: bfastfinal.fut helpers.fut mroc.fut recresid.fut
futhark pkg upgrade
futhark pkg sync
futhark pyopencl --library bfastfinal.fut
Loading

0 comments on commit a37c899

Please sign in to comment.