From d606f399db6cfa35f12cb64af16d000132feeb3a Mon Sep 17 00:00:00 2001 From: 12rambau Date: Wed, 30 Jun 2021 18:01:43 +0000 Subject: [PATCH 1/6] drop numpy for 2 items --- bfast/monitor/utils.py | 12 ++++-------- 1 file changed, 4 insertions(+), 8 deletions(-) diff --git a/bfast/monitor/utils.py b/bfast/monitor/utils.py index f69728c..81730d4 100644 --- a/bfast/monitor/utils.py +++ b/bfast/monitor/utils.py @@ -384,15 +384,12 @@ __critval_h = np.array([0.25, 0.5, 1]) __critval_period = np.arange(2, 12, 2) __critval_level = np.arange(0.95, 0.999, 0.001) -__critval_mr = np.array(["max", "range"]) +__critval_mr = ["max", "range"] def check(h, period, level, mr): if not h in __critval_h: raise ValueError("h can only be one of", __critval_h) - # if not period in __critval_period: - # raise ValueError("period can only be one of", __critval_period) - if not level in __critval_level: raise ValueError("level can only be one of", __critval_level) @@ -406,13 +403,12 @@ def get_critval(h, period, level, mr): index = np.zeros(4, dtype=np.int) # Get index into table from arguments - index[0] = np.where(mr == __critval_mr)[0][0] + index[0] = next(i for i, v in enumerate(__critval_mr) if v == mr) index[1] = np.where(level == __critval_level)[0][0] - # index[2] = np.where(period == __critval_period)[0][0] - # print((np.abs(__critval_period - period)).argmin()) index[2] = (np.abs(__critval_period - period)).argmin() index[3] = np.where(h == __critval_h)[0][0] - # For historical reasons, the critvals are scaled by sqrt(2) + + # For legacy reasons, the critvals are scaled by sqrt(2) return __critvals[tuple(index)] * np.sqrt(2) def _find_index_date(dates, t): From 13d984c7cb93f2c8773c26d8c9a499e1e46a6fc5 Mon Sep 17 00:00:00 2001 From: 12rambau Date: Wed, 30 Jun 2021 18:14:54 +0000 Subject: [PATCH 2/6] drop sklearn usage --- bfast/monitor/python/base.py | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/bfast/monitor/python/base.py b/bfast/monitor/python/base.py index ba496d7..17e87d4 100644 --- a/bfast/monitor/python/base.py +++ b/bfast/monitor/python/base.py @@ -10,7 +10,6 @@ import numpy as np np.warnings.filterwarnings('ignore') np.set_printoptions(suppress=True) -from sklearn import linear_model from bfast.base import BFASTMonitorBase from bfast.monitor.utils import compute_end_history, compute_lam, map_indices @@ -225,14 +224,14 @@ def fit_single(self, y): 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 - model = linear_model.LinearRegression(fit_intercept=False) - model.fit(X_nn_h.T, y_nn_h) + coef = np.linalg.pinv(X_nn_h@X_nn_h.T)@X_nn_h@y_nn_h + + if self.verbose > 1: - column_names = np.array(["Intercept", - "trend", + column_names = np.array(["trend", "harmonsin1", "harmoncos1", "harmonsin2", @@ -245,10 +244,10 @@ def fit_single(self, y): indxs = np.array([0, 2, 4, 6, 1, 3, 5]) # print(column_names[indxs]) print(column_names[indxs]) - print(model.coef_[indxs]) + print(coef[indxs]) # get predictions for all non-nan points - y_pred = model.predict(X_nn.T) + y_pred = X_nn.T@coef y_error = y_nn - y_pred # (2) evaluate model on monitoring period mosum_nn process From 69d21b0b5070993d9a87ed36e7abe7586c736d20 Mon Sep 17 00:00:00 2001 From: 12rambau Date: Thu, 1 Jul 2021 13:58:34 +0000 Subject: [PATCH 3/6] use nonzero instead of where https://numpy.org/doc/stable/reference/generated/numpy.where.html --- bfast/monitor/python/base.py | 2 +- bfast/monitor/utils.py | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/bfast/monitor/python/base.py b/bfast/monitor/python/base.py index 8cfe49b..a7b1b76 100644 --- a/bfast/monitor/python/base.py +++ b/bfast/monitor/python/base.py @@ -277,7 +277,7 @@ def fit_single(self, y): print("bounds", bounds) breaks = np.abs(mosum) > bounds - first_break = np.where(breaks)[0] + first_break = np.nonzero(breaks)[0] if first_break.shape[0] > 0: first_break = first_break[0] diff --git a/bfast/monitor/utils.py b/bfast/monitor/utils.py index 20853cf..bde3692 100644 --- a/bfast/monitor/utils.py +++ b/bfast/monitor/utils.py @@ -404,9 +404,9 @@ def get_critval(h, period, level, mr): # Get index into table from arguments index[0] = next(i for i, v in enumerate(__critval_mr) if v == mr) - index[1] = np.where(level == __critval_level)[0][0] + index[1] = np.nonzero(level == __critval_level)[0][0] index[2] = (np.abs(__critval_period - period)).argmin() - index[3] = np.where(h == __critval_h)[0][0] + index[3] = np.nonzero(h == __critval_h)[0][0] # For legacy reasons, the critvals are scaled by sqrt(2) return __critvals[tuple(index)] * np.sqrt(2) From 4e18a2ac89eb6762d0a18fa9fc03d27931bb26f9 Mon Sep 17 00:00:00 2001 From: 12rambau Date: Thu, 1 Jul 2021 13:59:39 +0000 Subject: [PATCH 4/6] remove chunks --- bfast/monitor/python/base.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/bfast/monitor/python/base.py b/bfast/monitor/python/base.py index a7b1b76..47aae29 100644 --- a/bfast/monitor/python/base.py +++ b/bfast/monitor/python/base.py @@ -101,7 +101,7 @@ def __init__(self, self._timers = {} self.use_mp = use_mp - def fit(self, data, dates, n_chunks=None, nan_value=0): + def fit(self, data, dates, nan_value=0, **kwargs): """ Fits the models for the ndarray 'data' Parameters @@ -114,7 +114,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 29b7ee77f7d608fd1246d0c6109f8628a216947b Mon Sep 17 00:00:00 2001 From: 12rambau Date: Thu, 1 Jul 2021 14:50:02 +0000 Subject: [PATCH 5/6] use apply_along_axis instead of for loop --- bfast/monitor/python/base.py | 48 +++++++++++------------------------- 1 file changed, 15 insertions(+), 33 deletions(-) diff --git a/bfast/monitor/python/base.py b/bfast/monitor/python/base.py index 47aae29..279f70f 100644 --- a/bfast/monitor/python/base.py +++ b/bfast/monitor/python/base.py @@ -150,30 +150,14 @@ def fit(self, data, dates, nan_value=0, **kwargs): self.magnitudes = rval[:,:,2].astype(np.float32) self.valids = rval[:,:,3].astype(np.int32) else: - means_global = np.zeros((data.shape[1], data.shape[2]), dtype=np.float32) - magnitudes_global = np.zeros((data.shape[1], data.shape[2]), dtype=np.float32) - breaks_global = np.zeros((data.shape[1], data.shape[2]), dtype=np.int32) - valids_global = np.zeros((data.shape[1], data.shape[2]), dtype=np.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 = 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 - self.means = means_global - self.magnitudes = magnitudes_global - self.valids = valids_global + + rval = np.apply_along_axis(self.fit_single, 0, data) + + #print(rval.shape) + self.breaks = rval[0].astype(np.int32) + self.means = rval[1].astype(np.float32) + self.magnitudes = rval[2].astype(np.float32) + self.valids = rval[3].astype(np.int32) return self @@ -227,21 +211,17 @@ def fit_single(self, y): # (1) fit linear regression model for history period coef = np.linalg.pinv(X_nn_h@X_nn_h.T)@X_nn_h@y_nn_h - - if self.verbose > 1: - column_names = np.array(["trend", - "harmonsin1", + column_names = np.array(["harmonsin1", "harmoncos1", "harmonsin2", "harmoncos2", "harmonsin3", "harmoncos3"]) if self.trend: - indxs = np.array([0, 1, 3, 5, 7, 2, 4, 6]) + indxs = np.array([1, 3, 5, 7, 2, 4, 6]) else: - indxs = np.array([0, 2, 4, 6, 1, 3, 5]) - # print(column_names[indxs]) + indxs = np.array([2, 4, 6, 1, 3, 5]) print(column_names[indxs]) print(coef[indxs]) @@ -279,11 +259,13 @@ def fit_single(self, y): first_break = np.nonzero(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, magnitude, Ns + rval = np.array([first_break, mean.item(), magnitude.item(), Ns.item()]) + + return rval def get_timers(self): """ Returns runtime measurements for the From 320cdd36426cc8d3457e2950902e545e55ccbf7d Mon Sep 17 00:00:00 2001 From: 12rambau Date: Fri, 2 Jul 2021 16:07:24 +0000 Subject: [PATCH 6/6] always return np.array --- bfast/monitor/python/base.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/bfast/monitor/python/base.py b/bfast/monitor/python/base.py index 279f70f..fd598c0 100644 --- a/bfast/monitor/python/base.py +++ b/bfast/monitor/python/base.py @@ -193,7 +193,10 @@ 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 + + rval = np.array([brk, mean, magnitude, Ns]) + + return rval val_inds = val_inds[ns:] val_inds -= self.n