From bbdf81f09188dcdd79a8dcc16e7cd551ff5e504f Mon Sep 17 00:00:00 2001 From: Pierrick Rambaud Date: Wed, 30 Jun 2021 16:01:52 +0000 Subject: [PATCH 1/9] first cupy implementation --- bfast/models.py | 13 ++ bfast/monitor/__init__.py | 1 + bfast/monitor/cupy/__init__.py | 1 + bfast/monitor/cupy/base.py | 296 +++++++++++++++++++++++++++++++++ bfast/monitor/cupy_utils.py | 87 ++++++++++ bfast/monitor/python/base.py | 2 +- 6 files changed, 399 insertions(+), 1 deletion(-) create mode 100644 bfast/monitor/cupy/__init__.py create mode 100644 bfast/monitor/cupy/base.py create mode 100644 bfast/monitor/cupy_utils.py diff --git a/bfast/models.py b/bfast/models.py index 98520bd..4c2b29a 100644 --- a/bfast/models.py +++ b/bfast/models.py @@ -1,5 +1,6 @@ from bfast.monitor import BFASTMonitorPython from bfast.monitor import BFASTMonitorOpenCL +from bfast.monitor import BFASTMonitorCuPy class BFASTMonitor(): @@ -151,6 +152,18 @@ def fit(self, data, dates, n_chunks=None, nan_value=0): period=self.period, verbose=self.verbose, ) + + elif self.backend == 'cupy': + self._model = BFASTMonitorCuPy( + start_monitor=self.start_monitor, + freq=self.freq, + k=self.k, + hfrac=self.hfrac, + trend=self.trend, + level=self.level, + period=self.period, + verbose=self.verbose, + ) elif self.backend == 'python-mp': self._model = BFASTMonitorPython( diff --git a/bfast/monitor/__init__.py b/bfast/monitor/__init__.py index 736de20..5f0c6fe 100644 --- a/bfast/monitor/__init__.py +++ b/bfast/monitor/__init__.py @@ -1,2 +1,3 @@ from .opencl import BFASTMonitorOpenCL from .python import BFASTMonitorPython +from .cupy import BFASTMonitorCuPy diff --git a/bfast/monitor/cupy/__init__.py b/bfast/monitor/cupy/__init__.py new file mode 100644 index 0000000..f0ac7c4 --- /dev/null +++ b/bfast/monitor/cupy/__init__.py @@ -0,0 +1 @@ +from .base import BFASTMonitorCuPy \ No newline at end of file diff --git a/bfast/monitor/cupy/base.py b/bfast/monitor/cupy/base.py new file mode 100644 index 0000000..18c145d --- /dev/null +++ b/bfast/monitor/cupy/base.py @@ -0,0 +1,296 @@ +''' +Created on June 15, 2021 + +@author: Pierrick Rambaud +''' + +import multiprocessing as mp +from functools import partial + +import cupy as cp +#cp.warnings.filterwarnings('ignore') +#cp.set_printoptions(suppress=True) +from sklearn import linear_model + +from bfast.base import BFASTMonitorBase +from bfast.monitor.cupy_utils import compute_end_history, compute_lam, map_indices + + +class BFASTMonitorCuPy(BFASTMonitorBase): + """ BFAST Monitor implementation based on Python and Numpy. The + interface follows the one of the corresponding R package + (see https://cran.r-project.org/web/packages/bfast) + + def __init__(self, + start_monitor, + freq=365, + k=3, + hfrac=0.25, + trend=True, + level=0.05, + detailed_results=False, + old_version=False, + verbose=0, + platform_id=0, + device_id=0 + ): + + Parameters + ---------- + + start_monitor : datetime object + A datetime object specifying the start of + the monitoring phase. + + freq : int, default 365 + The frequency for the seasonal model. + + k : int, default 3 + The number of harmonic terms. + + hfrac : float, default 0.25 + Float in the interval (0,1) specifying the + bandwidth relative to the sample size in + the MOSUM/ME monitoring processes. + + trend : bool, default True + Whether a tend offset term shall be used or not + + level : float, default 0.05 + Significance level of the monitoring (and ROC, + if selected) procedure, i.e., probability of + type I error. + + period : int, default 10 + Maximum time (relative to the history period) + that will be monitored. + + verbose : int, optional (default=0) + The verbosity level (0=no output, 1=output) + + use_mp : bool, default False + Determines whether to use the (very primitive) Python + multiprocessing or not. Enable for a speedup + + Examples + -------- + + Notes + ----- + + """ + def __init__(self, + start_monitor, + freq=365, + k=3, + hfrac=0.25, + trend=True, + level=0.05, + period=10, + verbose=0, + use_mp=False + ): + + super().__init__(start_monitor, + freq, + k=k, + hfrac=hfrac, + trend=trend, + level=level, + period=period, + verbose=verbose) + + self._timers = {} + self.use_mp = use_mp + + def fit(self, data, dates, n_chunks=None, nan_value=0): + """ Fits the models for the ndarray 'data' + + Parameters + ---------- + data: ndarray of shape (N, W, H), + where N is the number of time + series points per pixel and W + and H the width and the height + of the image, respectively. + dates : list of datetime objects + Specifies the dates of the elements + in data indexed by the first axis + n_chunks : int or None, default None + nan_value : int, default 0 + Specified the NaN value used in + the array data + + Returns + ------- + self : instance of BFASTMonitor + The object itself. + """ + data_ints = cp.copy(data) + data = cp.copy(data_ints).astype(cp.float32) + + # set NaN values + data[data_ints==nan_value] = cp.nan + + self.n = compute_end_history(dates, self.start_monitor) + + # create (complete) seasonal matrix ("patterns" as columns here!) + self.mapped_indices = map_indices(dates).astype(cp.int32) + self.X = self._create_data_matrix(self.mapped_indices) + + # period = data.shape[0] / cp.float(self.n) + self.lam = compute_lam(data.shape[0], self.hfrac, self.level, self.period) + + means_global = cp.zeros((data.shape[1], data.shape[2]), dtype=cp.float32) + magnitudes_global = cp.zeros((data.shape[1], data.shape[2]), dtype=cp.float32) + breaks_global = cp.zeros((data.shape[1], data.shape[2]), dtype=cp.int32) + valids_global = cp.zeros((data.shape[1], data.shape[2]), dtype=cp.int32) + + for i in range(data.shape[1]): + if self.verbose > 0: + print("Processing row {}".format(i)) + + for j in range(data.shape[2]): + y = cp.array(data[:,i,j]) + (pix_break, + pix_mean, + pix_magnitude, + pix_num_valid) = self.fit_single(y) + breaks_global[i,j] = pix_break + means_global[i,j] = pix_mean + magnitudes_global[i,j] = pix_magnitude + valids_global[i,j] = pix_num_valid + + self.breaks = breaks_global.get() + self.means = means_global.get() + self.magnitudes = magnitudes_global.get() + self.valids = valids_global.get() + + return self + + def fit_single(self, y): + """ Fits the BFAST model for the 1D array y. + + Parameters + ---------- + y : array + 1d array of length N + + Returns + ------- + self : instance of BFASTCPU + The object itself + """ + N = y.shape[0] + + # compute nan mappings + nans = cp.isnan(y) + num_nans = cp.cumsum(nans) + val_inds = cp.array(range(N))[~nans] + + # compute new limits (in data NOT containing missing values) + # ns = n - num_nans[self.n] + ns = self.n - num_nans[self.n - 1] + h = cp.int(float(ns) * self.hfrac) + Ns = N - num_nans[N - 1] + + if ns <= 5 or Ns - ns <= 5: + brk = -2 + mean = 0.0 + magnitude = 0.0 + if self.verbose > 1: + print("WARNING: Not enough observations: ns={ns}, Ns={Ns}".format(ns=ns, Ns=Ns)) + return brk, mean, magnitude, Ns + + val_inds = val_inds[ns:] + val_inds -= self.n + + # remove nan values from patterns+targets + X_nn = self.X[:, ~nans] + y_nn = y[~nans] + + # split data into history and monitoring phases + X_nn_h = X_nn[:, :ns] + X_nn_m = X_nn[:, ns:] + y_nn_h = y_nn[:ns] + y_nn_m = y_nn[ns:] + + # (1) fit linear regression model for history period + coef = cp.linalg.pinv(X_nn_h@X_nn_h.T)@X_nn_h@y_nn_h + + # get predictions for all non-nan points + y_pred = X_nn.T@coef + + y_error = y_nn - y_pred + + # (2) evaluate model on monitoring period mosum_nn process + err_cs = cp.cumsum(y_error[ns - h:Ns + 1]) + mosum_nn = err_cs[h:] - err_cs[:-h] + + sigma = cp.sqrt(cp.sum(y_error[:ns] ** 2) / (ns - (2 + 2 * self.k))) + mosum_nn = 1.0 / (sigma * cp.sqrt(ns)) * mosum_nn + + mosum = cp.full(N - self.n, cp.nan) + mosum[val_inds[:Ns - ns]] = mosum_nn + if self.verbose: + print("MOSUM process", mosum_nn.shape) + + # compute mean + mean = cp.mean(mosum_nn) + + # compute magnitude + magnitude = cp.median(y_error[ns:]) + + # boundary and breaks + a = self.mapped_indices[self.n:] / self.mapped_indices[self.n - 1].astype(cp.float) + bounds = self.lam * cp.sqrt(self._log_plus(a)) + + if self.verbose: + print("lambda", self.lam) + print("bounds", bounds) + + breaks = cp.abs(mosum) > bounds + first_break = cp.where(breaks)[0] + + if first_break.shape[0] > 0: + first_break = first_break[0] + else: + first_break = -1 + + return first_break, mean.get(), magnitude.get(), Ns + + def get_timers(self): + """ Returns runtime measurements for the + different phases of the fitting process. + + Returns + ------- + dict : An array containing the runtimes + for the different phases. + """ + return self._timers + + def _create_data_matrix(self, mapped_indices): + + # cast to cp array + mapped_indices = cp.array(mapped_indices) + N = mapped_indices.shape[0] + temp = 2 * cp.pi * mapped_indices / cp.float(self.freq) + + if self.trend: + X = cp.vstack((cp.ones(N), mapped_indices)) + else: + X = cp.ones(N) + + for j in cp.arange(1, self.k + 1): + X = cp.vstack((X, cp.sin(j * temp))) + X = cp.vstack((X, cp.cos(j * temp))) + + return X + + def _log_plus(self, a): + retval = cp.ones(a.shape, dtype=cp.float) + fl = a > cp.e + retval[fl] = cp.log(a[fl]) + + return retval diff --git a/bfast/monitor/cupy_utils.py b/bfast/monitor/cupy_utils.py new file mode 100644 index 0000000..87d72f2 --- /dev/null +++ b/bfast/monitor/cupy_utils.py @@ -0,0 +1,87 @@ +# copy the utils functionalities into cupy objects +from datetime import datetime + +import cupy as cp +import numpy as np +import pandas + +from . import utils + + +__critvals = cp.array(utils.__critvals) +__critval_h = cp.array(utils.__critval_h) +__critval_period = cp.array(utils.__critval_period) +__critval_level = cp.array(utils.__critval_level) +__critval_mr = utils.__critval_mr.tolist() + +check = utils.check + +def get_critval(h, period, level, mr): + + # Sanity check + check(h, period, level, mr) + + index = cp.zeros(4, dtype=cp.int) + + # Get index into table from arguments + index[0] = next(i for i, v in enumerate(__critval_mr) if v == mr) + index[1] = cp.where(level == __critval_level)[0][0] + index[2] = (cp.abs(__critval_period - period)).argmin() + index[3] = cp.where(h == __critval_h)[0][0] + + # For legacy reasons, the critvals are scaled by sqrt(2) + return __critvals[tuple(index)] * cp.sqrt(2) + +_find_index_date = utils._find_index_date + +def crop_data_dates(data, dates, start, end): + """ Crops the input data and the associated + dates w.r.t. the provided start and end + datetime object. + + Parameters + ---------- + data: ndarray of shape (N, W, H) + Here, N is the number of time + series points per pixel and W + and H are the width and the height + of the image, respectively. + dates : list of datetime objects + Specifies the dates of the elements + in data indexed by the first axis + n_chunks : int or None, default None + start : datetime + The start datetime object + end : datetime + The end datetime object + + Returns + ------- + Returns: data, dates + The cropped data array and the + cropped list. Only those images + and dates that are with the start/end + period are contained in the returned + objects. + """ + + start_idx = _find_index_date(dates, start) + end_idx = _find_index_date(dates, end) + + data_cropped = data[start_idx:end_idx, :, :] + dates_cropped = list(cp.array(dates)[start_idx:end_idx]) + + return data_cropped, dates_cropped + + +def compute_lam(N, hfrac, level, period): + + check(hfrac, period, 1 - level, "max") + + return get_critval(hfrac, period, 1 - level, "max") + +compute_end_history = utils.compute_end_history + +def map_indices(dates): + + return cp.array(utils.map_indices(dates)) diff --git a/bfast/monitor/python/base.py b/bfast/monitor/python/base.py index ba496d7..5a6d593 100644 --- a/bfast/monitor/python/base.py +++ b/bfast/monitor/python/base.py @@ -258,7 +258,7 @@ def fit_single(self, y): sigma = np.sqrt(np.sum(y_error[:ns] ** 2) / (ns - (2 + 2 * self.k))) mosum_nn = 1.0 / (sigma * np.sqrt(ns)) * mosum_nn - mosum = np.repeat(np.nan, N - self.n) + mosum = np.full(N - self.n, np.nan) mosum[val_inds[:Ns - ns]] = mosum_nn if self.verbose: print("MOSUM process", mosum_nn.shape) From a872a2425c9a60194508e8fa0567488a9202cdf5 Mon Sep 17 00:00:00 2001 From: 12rambau Date: Thu, 1 Jul 2021 06:09:36 +0000 Subject: [PATCH 2/9] select the GPU device --- bfast/monitor/cupy/base.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/bfast/monitor/cupy/base.py b/bfast/monitor/cupy/base.py index 18c145d..4e397f6 100644 --- a/bfast/monitor/cupy/base.py +++ b/bfast/monitor/cupy/base.py @@ -88,7 +88,7 @@ def __init__(self, level=0.05, period=10, verbose=0, - use_mp=False + device_id=0 ): super().__init__(start_monitor, @@ -101,7 +101,9 @@ def __init__(self, verbose=verbose) self._timers = {} - self.use_mp = use_mp + + # use the specified cuda device (default to 0) + cp.cuda.device.Device(device_id).use() def fit(self, data, dates, n_chunks=None, nan_value=0): """ Fits the models for the ndarray 'data' @@ -126,6 +128,7 @@ def fit(self, data, dates, n_chunks=None, nan_value=0): self : instance of BFASTMonitor The object itself. """ + data_ints = cp.copy(data) data = cp.copy(data_ints).astype(cp.float32) @@ -140,7 +143,7 @@ def fit(self, data, dates, n_chunks=None, nan_value=0): # period = data.shape[0] / cp.float(self.n) self.lam = compute_lam(data.shape[0], self.hfrac, self.level, self.period) - + means_global = cp.zeros((data.shape[1], data.shape[2]), dtype=cp.float32) magnitudes_global = cp.zeros((data.shape[1], data.shape[2]), dtype=cp.float32) breaks_global = cp.zeros((data.shape[1], data.shape[2]), dtype=cp.int32) From f01cfc9836dfc9dec5ea4c32fb44efc2e61b1101 Mon Sep 17 00:00:00 2001 From: 12rambau Date: Thu, 1 Jul 2021 06:21:55 +0000 Subject: [PATCH 3/9] drop n_chunks --- bfast/monitor/cupy/base.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/bfast/monitor/cupy/base.py b/bfast/monitor/cupy/base.py index 4e397f6..e96d2b5 100644 --- a/bfast/monitor/cupy/base.py +++ b/bfast/monitor/cupy/base.py @@ -105,7 +105,7 @@ def __init__(self, # use the specified cuda device (default to 0) cp.cuda.device.Device(device_id).use() - def fit(self, data, dates, n_chunks=None, nan_value=0): + def fit(self, data, dates, nan_value=0): """ Fits the models for the ndarray 'data' Parameters @@ -118,7 +118,6 @@ def fit(self, data, dates, n_chunks=None, nan_value=0): dates : list of datetime objects Specifies the dates of the elements in data indexed by the first axis - n_chunks : int or None, default None nan_value : int, default 0 Specified the NaN value used in the array data From 0e45fb97c259c09ace6ac3a36ae3f998d0d89fdb Mon Sep 17 00:00:00 2001 From: 12rambau Date: Thu, 1 Jul 2021 13:54:58 +0000 Subject: [PATCH 4/9] test along_axis --- bfast/monitor/cupy/base.py | 72 +++++++++++++++++++++--------------- bfast/monitor/python/base.py | 9 ++++- 2 files changed, 51 insertions(+), 30 deletions(-) diff --git a/bfast/monitor/cupy/base.py b/bfast/monitor/cupy/base.py index e96d2b5..40310cd 100644 --- a/bfast/monitor/cupy/base.py +++ b/bfast/monitor/cupy/base.py @@ -105,7 +105,7 @@ def __init__(self, # use the specified cuda device (default to 0) cp.cuda.device.Device(device_id).use() - def fit(self, data, dates, nan_value=0): + def fit(self, data, dates, nan_value=0, **kwargs): """ Fits the models for the ndarray 'data' Parameters @@ -129,7 +129,7 @@ def fit(self, data, dates, nan_value=0): """ data_ints = cp.copy(data) - data = cp.copy(data_ints).astype(cp.float32) + data = cp.array(cp.copy(data_ints)).astype(cp.float32) # set NaN values data[data_ints==nan_value] = cp.nan @@ -143,30 +143,40 @@ def fit(self, data, dates, nan_value=0): # period = data.shape[0] / cp.float(self.n) self.lam = compute_lam(data.shape[0], self.hfrac, self.level, self.period) - means_global = cp.zeros((data.shape[1], data.shape[2]), dtype=cp.float32) - magnitudes_global = cp.zeros((data.shape[1], data.shape[2]), dtype=cp.float32) - breaks_global = cp.zeros((data.shape[1], data.shape[2]), dtype=cp.int32) - valids_global = cp.zeros((data.shape[1], data.shape[2]), dtype=cp.int32) - - for i in range(data.shape[1]): - if self.verbose > 0: - print("Processing row {}".format(i)) - - for j in range(data.shape[2]): - y = cp.array(data[:,i,j]) - (pix_break, - pix_mean, - pix_magnitude, - pix_num_valid) = self.fit_single(y) - breaks_global[i,j] = pix_break - means_global[i,j] = pix_mean - magnitudes_global[i,j] = pix_magnitude - valids_global[i,j] = pix_num_valid - - self.breaks = breaks_global.get() - self.means = means_global.get() - self.magnitudes = magnitudes_global.get() - self.valids = valids_global.get() + #means_global = cp.zeros((data.shape[1], data.shape[2]), dtype=cp.float32) + #magnitudes_global = cp.zeros((data.shape[1], data.shape[2]), dtype=cp.float32) + #breaks_global = cp.zeros((data.shape[1], data.shape[2]), dtype=cp.int32) + #valids_global = cp.zeros((data.shape[1], data.shape[2]), dtype=cp.int32) + print(data.shape) + + results = cp.apply_along_axis(self.fit_single, axis=0, arr=data) + print(results) + print(results.shape) + + self.breaks = results[:,:,0].get() + self.means = results[:,:,1].get() + self.magnitudes = results[:,:,2].get() + self.valids = results[:,:,3].get() + + #for i in range(data.shape[1]): + # if self.verbose > 0: + # print("Processing row {}".format(i)) +# + # for j in range(data.shape[2]): + # y = cp.array(data[:,i,j]) + # (pix_break, + # pix_mean, + # pix_magnitude, + # pix_num_valid) = self.fit_single(y) + # breaks_global[i,j] = pix_break + # means_global[i,j] = pix_mean + # magnitudes_global[i,j] = pix_magnitude + # valids_global[i,j] = pix_num_valid +# + #self.breaks = breaks_global.get() + #self.means = means_global.get() + #self.magnitudes = magnitudes_global.get() + #self.valids = valids_global.get() return self @@ -183,6 +193,8 @@ def fit_single(self, y): self : instance of BFASTCPU The object itself """ + print(y) + print('toto') N = y.shape[0] # compute nan mappings @@ -255,11 +267,13 @@ def fit_single(self, y): first_break = cp.where(breaks)[0] if first_break.shape[0] > 0: - first_break = first_break[0] + first_break = first_break[0].item() else: first_break = -1 - - return first_break, mean.get(), magnitude.get(), Ns + + res = cp.array([first_break, mean.item(), magnitude.item(), Ns.item()]) + + print(res) def get_timers(self): """ Returns runtime measurements for the diff --git a/bfast/monitor/python/base.py b/bfast/monitor/python/base.py index 5a6d593..729e3e8 100644 --- a/bfast/monitor/python/base.py +++ b/bfast/monitor/python/base.py @@ -163,15 +163,22 @@ def fit(self, data, dates, n_chunks=None, nan_value=0): for j in range(data.shape[2]): y = data[:,i,j] + + print(y.shape) + (pix_break, pix_mean, pix_magnitude, pix_num_valid) = self.fit_single(y) + print(self.fit_single(y)) + breaks_global[i,j] = pix_break means_global[i,j] = pix_mean magnitudes_global[i,j] = pix_magnitude valids_global[i,j] = pix_num_valid - + + raise Exception('stop') + self.breaks = breaks_global self.means = means_global self.magnitudes = magnitudes_global From cbb3372979956acc5605a97c485d166316af3671 Mon Sep 17 00:00:00 2001 From: 12rambau Date: Fri, 2 Jul 2021 14:08:13 +0000 Subject: [PATCH 5/9] prevent cupy to load if no GPU is availabel --- bfast/models.py | 7 ++----- bfast/monitor/__init__.py | 12 +++++++++++- bfast/monitor/cupy/base.py | 37 ++++++------------------------------- 3 files changed, 19 insertions(+), 37 deletions(-) diff --git a/bfast/models.py b/bfast/models.py index 4c2b29a..b671d2e 100644 --- a/bfast/models.py +++ b/bfast/models.py @@ -1,7 +1,4 @@ -from bfast.monitor import BFASTMonitorPython -from bfast.monitor import BFASTMonitorOpenCL -from bfast.monitor import BFASTMonitorCuPy - +from bfast.monitor import * class BFASTMonitor(): """ @@ -153,7 +150,7 @@ def fit(self, data, dates, n_chunks=None, nan_value=0): verbose=self.verbose, ) - elif self.backend == 'cupy': + elif self.backend == 'cupy' and gpu_available: self._model = BFASTMonitorCuPy( start_monitor=self.start_monitor, freq=self.freq, diff --git a/bfast/monitor/__init__.py b/bfast/monitor/__init__.py index 5f0c6fe..a816f42 100644 --- a/bfast/monitor/__init__.py +++ b/bfast/monitor/__init__.py @@ -1,3 +1,13 @@ from .opencl import BFASTMonitorOpenCL from .python import BFASTMonitorPython -from .cupy import BFASTMonitorCuPy + +# only import cupy related bfast if a gpu is available +try: + import cupy + cupy.cuda.Device() + gpu_available = True +except: + gpu_available = False + +if gpu_available: + from .cupy import BFASTMonitorCuPy diff --git a/bfast/monitor/cupy/base.py b/bfast/monitor/cupy/base.py index 40310cd..110be86 100644 --- a/bfast/monitor/cupy/base.py +++ b/bfast/monitor/cupy/base.py @@ -143,40 +143,16 @@ def fit(self, data, dates, nan_value=0, **kwargs): # period = data.shape[0] / cp.float(self.n) self.lam = compute_lam(data.shape[0], self.hfrac, self.level, self.period) - #means_global = cp.zeros((data.shape[1], data.shape[2]), dtype=cp.float32) - #magnitudes_global = cp.zeros((data.shape[1], data.shape[2]), dtype=cp.float32) - #breaks_global = cp.zeros((data.shape[1], data.shape[2]), dtype=cp.int32) - #valids_global = cp.zeros((data.shape[1], data.shape[2]), dtype=cp.int32) print(data.shape) results = cp.apply_along_axis(self.fit_single, axis=0, arr=data) print(results) print(results.shape) - self.breaks = results[:,:,0].get() - self.means = results[:,:,1].get() - self.magnitudes = results[:,:,2].get() - self.valids = results[:,:,3].get() - - #for i in range(data.shape[1]): - # if self.verbose > 0: - # print("Processing row {}".format(i)) -# - # for j in range(data.shape[2]): - # y = cp.array(data[:,i,j]) - # (pix_break, - # pix_mean, - # pix_magnitude, - # pix_num_valid) = self.fit_single(y) - # breaks_global[i,j] = pix_break - # means_global[i,j] = pix_mean - # magnitudes_global[i,j] = pix_magnitude - # valids_global[i,j] = pix_num_valid -# - #self.breaks = breaks_global.get() - #self.means = means_global.get() - #self.magnitudes = magnitudes_global.get() - #self.valids = valids_global.get() + self.breaks = results[0].get() + self.means = results[1].get() + self.magnitudes = results[2].get() + self.valids = results[3].get() return self @@ -193,8 +169,7 @@ def fit_single(self, y): self : instance of BFASTCPU The object itself """ - print(y) - print('toto') + N = y.shape[0] # compute nan mappings @@ -273,7 +248,7 @@ def fit_single(self, y): res = cp.array([first_break, mean.item(), magnitude.item(), Ns.item()]) - print(res) + return res def get_timers(self): """ Returns runtime measurements for the From 342391c9e887305d680dfa1bfc8e06e47fad0d69 Mon Sep 17 00:00:00 2001 From: 12rambau Date: Fri, 2 Jul 2021 17:05:41 +0000 Subject: [PATCH 6/9] ensure consistensy of outputs --- bfast/monitor/cupy/base.py | 11 +++-------- 1 file changed, 3 insertions(+), 8 deletions(-) diff --git a/bfast/monitor/cupy/base.py b/bfast/monitor/cupy/base.py index 110be86..580dbf2 100644 --- a/bfast/monitor/cupy/base.py +++ b/bfast/monitor/cupy/base.py @@ -142,12 +142,8 @@ def fit(self, data, dates, nan_value=0, **kwargs): # period = data.shape[0] / cp.float(self.n) self.lam = compute_lam(data.shape[0], self.hfrac, self.level, self.period) - - print(data.shape) results = cp.apply_along_axis(self.fit_single, axis=0, arr=data) - print(results) - print(results.shape) self.breaks = results[0].get() self.means = results[1].get() @@ -189,7 +185,8 @@ def fit_single(self, y): magnitude = 0.0 if self.verbose > 1: print("WARNING: Not enough observations: ns={ns}, Ns={Ns}".format(ns=ns, Ns=Ns)) - return brk, mean, magnitude, Ns + + return cp.array([brk, mean, magnitude, Ns], dtype=cp.float) val_inds = val_inds[ns:] val_inds -= self.n @@ -246,9 +243,7 @@ def fit_single(self, y): else: first_break = -1 - res = cp.array([first_break, mean.item(), magnitude.item(), Ns.item()]) - - return res + return cp.array([first_break, mean.item(), magnitude.item(), Ns.item()], dtype=cp.float) def get_timers(self): """ Returns runtime measurements for the From d384b3b9ada17b47cd7fecde0b82127111dcfd98 Mon Sep 17 00:00:00 2001 From: 12rambau Date: Wed, 7 Jul 2021 13:25:59 +0000 Subject: [PATCH 7/9] drop cupy_utils --- bfast/monitor/cupy/base.py | 6 +-- bfast/monitor/cupy_utils.py | 87 ------------------------------------- 2 files changed, 3 insertions(+), 90 deletions(-) delete mode 100644 bfast/monitor/cupy_utils.py diff --git a/bfast/monitor/cupy/base.py b/bfast/monitor/cupy/base.py index 580dbf2..8f593f1 100644 --- a/bfast/monitor/cupy/base.py +++ b/bfast/monitor/cupy/base.py @@ -13,7 +13,7 @@ from sklearn import linear_model from bfast.base import BFASTMonitorBase -from bfast.monitor.cupy_utils import compute_end_history, compute_lam, map_indices +from bfast.monitor.utils import compute_end_history, compute_lam, map_indices class BFASTMonitorCuPy(BFASTMonitorBase): @@ -137,11 +137,11 @@ def fit(self, data, dates, nan_value=0, **kwargs): self.n = compute_end_history(dates, self.start_monitor) # create (complete) seasonal matrix ("patterns" as columns here!) - self.mapped_indices = map_indices(dates).astype(cp.int32) + self.mapped_indices = np.array(map_indices(dates)).astype(cp.int32) self.X = self._create_data_matrix(self.mapped_indices) # period = data.shape[0] / cp.float(self.n) - self.lam = compute_lam(data.shape[0], self.hfrac, self.level, self.period) + self.lam = cp.array(compute_lam(data.shape[0], self.hfrac, self.level, self.period)) results = cp.apply_along_axis(self.fit_single, axis=0, arr=data) diff --git a/bfast/monitor/cupy_utils.py b/bfast/monitor/cupy_utils.py deleted file mode 100644 index 87d72f2..0000000 --- a/bfast/monitor/cupy_utils.py +++ /dev/null @@ -1,87 +0,0 @@ -# copy the utils functionalities into cupy objects -from datetime import datetime - -import cupy as cp -import numpy as np -import pandas - -from . import utils - - -__critvals = cp.array(utils.__critvals) -__critval_h = cp.array(utils.__critval_h) -__critval_period = cp.array(utils.__critval_period) -__critval_level = cp.array(utils.__critval_level) -__critval_mr = utils.__critval_mr.tolist() - -check = utils.check - -def get_critval(h, period, level, mr): - - # Sanity check - check(h, period, level, mr) - - index = cp.zeros(4, dtype=cp.int) - - # Get index into table from arguments - index[0] = next(i for i, v in enumerate(__critval_mr) if v == mr) - index[1] = cp.where(level == __critval_level)[0][0] - index[2] = (cp.abs(__critval_period - period)).argmin() - index[3] = cp.where(h == __critval_h)[0][0] - - # For legacy reasons, the critvals are scaled by sqrt(2) - return __critvals[tuple(index)] * cp.sqrt(2) - -_find_index_date = utils._find_index_date - -def crop_data_dates(data, dates, start, end): - """ Crops the input data and the associated - dates w.r.t. the provided start and end - datetime object. - - Parameters - ---------- - data: ndarray of shape (N, W, H) - Here, N is the number of time - series points per pixel and W - and H are the width and the height - of the image, respectively. - dates : list of datetime objects - Specifies the dates of the elements - in data indexed by the first axis - n_chunks : int or None, default None - start : datetime - The start datetime object - end : datetime - The end datetime object - - Returns - ------- - Returns: data, dates - The cropped data array and the - cropped list. Only those images - and dates that are with the start/end - period are contained in the returned - objects. - """ - - start_idx = _find_index_date(dates, start) - end_idx = _find_index_date(dates, end) - - data_cropped = data[start_idx:end_idx, :, :] - dates_cropped = list(cp.array(dates)[start_idx:end_idx]) - - return data_cropped, dates_cropped - - -def compute_lam(N, hfrac, level, period): - - check(hfrac, period, 1 - level, "max") - - return get_critval(hfrac, period, 1 - level, "max") - -compute_end_history = utils.compute_end_history - -def map_indices(dates): - - return cp.array(utils.map_indices(dates)) From 3011509a1ef8036386d7d1d08519fcdd0f47758e Mon Sep 17 00:00:00 2001 From: 12rambau Date: Wed, 7 Jul 2021 13:31:31 +0000 Subject: [PATCH 8/9] clean-up documentation --- bfast/models.py | 7 ++++--- bfast/monitor/cupy/base.py | 19 +++---------------- 2 files changed, 7 insertions(+), 19 deletions(-) diff --git a/bfast/models.py b/bfast/models.py index b671d2e..a68dd77 100644 --- a/bfast/models.py +++ b/bfast/models.py @@ -5,7 +5,8 @@ class BFASTMonitor(): BFASTMonitor implements the BFASTMonitor approach and provides two backends/implementations: - - A pure Python implementation (based on the Numpy package). + - A pure Numpy implementation (based on the Numpy package) + - A pure Cupy implementation suited for GPU devices - An optimized OpenCL implementation suited for massively-parallel devices. The interface follows the one of the corresponding R package, @@ -50,8 +51,8 @@ class BFASTMonitor(): Specifies the OpenCL platform id. device_id int, default 0 - Only relevant if backend='opencl'. - Specified the OpenCL device id. + Only relevant if backend='opencl' or 'cupy'. + Specified the GPU device id. detailed_results : bool, default False Only relevant if backend='opencl'. diff --git a/bfast/monitor/cupy/base.py b/bfast/monitor/cupy/base.py index 8f593f1..a4b39d7 100644 --- a/bfast/monitor/cupy/base.py +++ b/bfast/monitor/cupy/base.py @@ -8,8 +8,6 @@ from functools import partial import cupy as cp -#cp.warnings.filterwarnings('ignore') -#cp.set_printoptions(suppress=True) from sklearn import linear_model from bfast.base import BFASTMonitorBase @@ -28,10 +26,7 @@ def __init__(self, hfrac=0.25, trend=True, level=0.05, - detailed_results=False, - old_version=False, verbose=0, - platform_id=0, device_id=0 ): @@ -67,17 +62,9 @@ def __init__(self, verbose : int, optional (default=0) The verbosity level (0=no output, 1=output) - - use_mp : bool, default False - Determines whether to use the (very primitive) Python - multiprocessing or not. Enable for a speedup - - Examples - -------- - - Notes - ----- - + + device_id : int, optional (default=1) + The GPU device id number """ def __init__(self, start_monitor, From 19c4084952228f8349618738e291e37e7ac590e5 Mon Sep 17 00:00:00 2001 From: 12rambau Date: Wed, 7 Jul 2021 13:33:31 +0000 Subject: [PATCH 9/9] add cupy to requirements --- requirements.txt | 1 + setup.py | 1 + 2 files changed, 2 insertions(+) diff --git a/requirements.txt b/requirements.txt index ad7e054..89e7d92 100644 --- a/requirements.txt +++ b/requirements.txt @@ -8,3 +8,4 @@ wget==3.2 Sphinx==2.2.0 sphinx-bootstrap-theme==0.7.1 numpydoc==1.0.0 +cupy==9.2.0 diff --git a/setup.py b/setup.py index 2426520..1a4ff26 100644 --- a/setup.py +++ b/setup.py @@ -93,6 +93,7 @@ def setup_package(): 'scipy>=1.2.1', 'matplotlib>=2.2.2', 'wget>=3.2', + 'cupy>=9.2.0' ], classifiers=['Intended Audience :: Science/Research', 'Intended Audience :: Developers',