# # ABOUT # Artisan Main Canvas # LICENSE # This program or module is free software: you can redistribute it and/or # modify it under the terms of the GNU General Public License as published # by the Free Software Foundation, either version 2 of the License, or # version 3 of the License, or (at your option) any later version. It is # provided for educational purposes and is distributed in the hope that # it will be useful, but WITHOUT ANY WARRANTY; without even the implied # warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See # the GNU General Public License for more details. # AUTHOR # Marko Luther, 2023 # SPDX-License-Identifier: GPL-3.0-or-later # Extracted from Artisan v4.0.2, commit 1f678a8336c28c4dc0ac08b9a28eccbc361130e5. # Original copyright: Artisan contributors / Marko Luther. GPL v2 or later; # this local derivative selects GPL v3 or later. Full license ships alongside. # This file contains calculation functions only; no device or control functions. from __future__ import annotations import numpy import numpy.typing as npt import math,sys,warnings,logging from typing import cast,Literal,Sequence class QApplication: @staticmethod def translate(_context, text): return text _log=logging.getLogger(__name__) class Artisan402Core: def inputFilter(self, timex:list[float], tempx:list[float], time:float, temp:float, BT:bool = False) -> float: try: wrong_reading = 0 ######################### # a) detect duplicates: remove a reading if it is equal to the previous or if that is -1 to the one before that if self.dropDuplicates and ((len(tempx)>1 and tempx[-1] == -1 and abs(temp - tempx[-2]) <= self.dropDuplicatesLimit) or (len(tempx)>0 and abs(temp - tempx[-1]) <= self.dropDuplicatesLimit)): wrong_reading = 2 # replace by previous reading not by -1 ######################### # b) detect overflows if self.minmaxLimits and (temp < self.filterDropOut_tmin or temp > self.filterDropOut_tmax): wrong_reading = 1 ######################### # c) detect spikes (on BT only after CHARGE if autoChargeFlag=True not to have a conflict here) n = self.filterDropOut_spikeRoR_period dRoR_limit = self.filterDropOut_spikeRoR_dRoR_limit # the limit of additional RoR in temp/sec (4C for C / 7F for F) compared to previous readings if self.dropSpikes and ((not self.autoChargeFlag) or (not BT) or (self.timeindex[0] != -1 and (self.timeindex[0] + n) < len(timex))) and not wrong_reading and len(tempx) >= n: # no min/max overflow detected # check if RoR caused by actual measurement is way higher then the previous one # calc previous RoR (pRoR) taking the last n samples into account pdtemp = tempx[-1] - tempx[-n] pdtime = timex[-1] - timex[-n] if pdtime > 0: pRoR = abs(pdtemp/pdtime) dtemp = tempx[-1] - temp dtime = timex[-1] - time if dtime > 0: RoR = abs(dtemp/dtime) if RoR > (pRoR + dRoR_limit): wrong_reading = 2 ######################### # c) handle outliers if it could be detected if wrong_reading: if len(tempx) > 0 and tempx[-1] != -1: # repeate last correct reading if not done before in the last two fixes (min/max violation are always filtered) if len(tempx) == 1 or (len(tempx) > 3 and (tempx[-1] != tempx[-2] or tempx[-2] != tempx[-3])): return tempx[-1] if wrong_reading == 1: return -1 # no way to correct this return temp if wrong_reading == 1: return -1 # no way to correct this return temp # try to improve a previously corrected reading timex/temp[-1] based on the current reading time/temp (just in this case the actual reading is not a drop) if (self.minmaxLimits or self.dropSpikes or self.dropDuplicates): if len(tempx) > 2 and tempx[-1] == tempx[-2] == tempx[-3] and tempx[-1] != -1 and tempx[-1] != temp and temp!=-1: # previous reading was a drop and replaced by reading[-2] and same for the one before delta = (tempx[-3] - temp) / 3.0 tempx[-1] = tempx[-3] - 2*delta tempx[-2] = tempx[-3] - delta elif len(tempx) > 1 and tempx[-1] == tempx[-2] and tempx[-1] != -1 and tempx[-1] != temp and temp!=-1: # previous reading was a drop and replaced by reading[-2] tempx[-1] = (tempx[-2] + temp) / 2.0 return temp except Exception as e: # pylint: disable=broad-except _log.exception(e) _, _, exc_tb = sys.exc_info() self.adderror((QApplication.translate('Error Message','Exception:') + ' filterDropOuts() {0}').format(str(e)),getattr(exc_tb, 'tb_lineno', '?')) return temp def decay_average(self, tx_in:list[float], temp_in:Sequence[float|None], decay_weights:list[int]|None) -> float: if decay_weights is None or len(decay_weights)<2 or len(tx_in) != len(temp_in): if len(temp_in)>0 and temp_in[-1] is not None: return temp_in[-1] # pyrefly: ignore[bad-return] return -1 l = min(len(decay_weights),len(temp_in)) # take trail of length l and remove items where temp[i]=None to fulfil precond. of numpy.interp tx_org:list[float] = [] temp_trail:list[float] = [] for x, tp in zip(tx_in[-l:], temp_in[-l:], strict=True): # we only iterate over l-elements if tp is not None and tp != -1: tx_org.append(x) temp_trail.append(tp) if len(temp_trail) == 0: # no valid values return -1 l = len(temp_trail) # might be shorter than before # len(tx)=len(temp) here and it is guaranteed that len(tx_org)=len(temp_trail) = l d = self.delay / 1000. # we create a linearly spaced time array starting from the newest timestamp in sampling interval distance tx_lin = numpy.flip(numpy.arange(tx_org[-1],tx_org[-1]-l*d,-d), axis=0) # by construction, len(tx_lin)=len(tx_org)=l temp_trail_re = numpy.interp(tx_lin, tx_org, temp_trail) # resample data into that linear spaced time try: return float(numpy.average(temp_trail_re[-len(decay_weights):],axis=0,weights=decay_weights[-l:])) except Exception: # pylint: disable=broad-except # in case something goes very wrong we at least return the standard average over temp, this should always work as len(tx)=len(temp) return float(numpy.average(numpy.array(temp_trail))) def smooth(self, x:'npt.NDArray[numpy.floating]', y:'npt.NDArray[numpy.float64]', window_len:int = 15, window:str = 'hanning') -> 'npt.NDArray[numpy.floating]': try: if len(x) == len(y) and len(x) > 1: if window_len > 2: # smooth curves #s = numpy.r_[2*x[0]-y[window_len:1:-1],y,2*y[-1]-y[-1:-window_len:-1]] #s=numpy.r_[y[window_len-1:0:-1],y,y[-2:-window_len-1:-1]] #s = y s = numpy.r_[y[window_len-1:0:-1],y,y[-1:-window_len:-1]] if window == 'flat': #moving average w = numpy.ones(window_len,'d') else: w = eval('numpy.'+window+'(window_len)') # pylint: disable=eval-used try: ys = numpy.convolve(w/w.sum(), s, mode='valid') except Exception: # pylint: disable=broad-except return y hwl = int(window_len/2) res = ys[hwl:-hwl] if len(res)+1 == len(y) and len(res) > 0: try: return ys[hwl-1:-hwl] except Exception: # pylint: disable=broad-except return y elif len(res) != len(y): return y return res return y return y except Exception as ex: # pylint: disable=broad-except _log.exception(ex) _, _, exc_tb = sys.exc_info() self.adderror((QApplication.translate('Error Message','Exception:') + ' smooth() {0}').format(str(ex)),getattr(exc_tb, 'tb_lineno', '?')) return x def smooth_slice(self, a:'npt.NDArray[numpy.double]', b:'npt.NDArray[numpy.float64]', window_len:int = 7, window:str = 'hanning', decay_weights:list[int]|None = None, decay_smoothing:bool = False, re_sample:bool = True, back_sample:bool = True, a_lin:'npt.NDArray[numpy.double]|None' = None, delta:bool=False) -> 'npt.NDArray[numpy.double]': a_mod:npt.NDArray[numpy.floating] # 1. re-sample if re_sample: if a_lin is None or len(a_lin) != len(a): a_mod = cast(numpy.ndarray[tuple[Literal[1]]], numpy.linspace(a[0],a[-1],len(a))) else: a_mod = a_lin b = numpy.interp(a_mod, a, b) # resample data to linear spaced time else: a_mod = a res:npt.NDArray[numpy.floating] = b # just in case the precondition (self.filterDropOuts or window_len>2) does not hold # 2. filter spikes (only applied offline) if self.filterDropOuts and not self.flagon: try: # if self.flagon: # online_medfilt = LiveMedian(median_filter_factor) # b = numpy.array(list(map(online_medfilt, b))) bb = self.medfilt(b, (self.median_filter_factor_RoR if delta else self.median_filter_factor)) # scipyernative which performs equal, but produces larger artefacts at the borders and for intermediate NaN values for k>3 # from scipy.signal import medfilt as scipy_medfilt # b = scipy_medfilt(b,3) res = bb except Exception as e: # pylint: disable=broad-except _log.exception(e) res = b # 3. smooth data if window_len>2: if decay_smoothing: # decay smoothing decay_weights_internal:npt.NDArray[numpy.int_] if decay_weights is None: decay_weights_internal = numpy.arange(1,window_len+1) else: window_len = len(decay_weights) decay_weights_internal = numpy.array(decay_weights) # invariant: window_len = len(decay_weights_internal) if decay_weights_internal.sum() == 0: res = b else: result:list[float] = [] # ignore -1 readings in averaging and ensure a good ramp for i, v in enumerate(b): seq = b[max(0,i-window_len + 1):i+1] w = decay_weights_internal[max(0,window_len-len(seq)):] # preCond: len(decay_weights_internal)=window_len and len(seq) <= window_len; postCond: len(w)=len(seq) if len(w) == 0: # we don't average if there is are no weights (e.g. if the original seq did only contain -1 values and got empty) result.append(v) else: result.append(float(numpy.average(seq,axis=0,weights=w))) # works only if len(seq) = len(w) res = numpy.array(result) # postCond: len(res) = len(b) else: # optimal smoothing (the default) win_len = max(0,window_len) if win_len != 1: # at the lowest level we turn smoothing completely off res = self.smooth(a_mod, b, win_len, window) # pyright:ignore[reportUnknownArgumentType] else: res = b # 4. sample back if re_sample and back_sample: res = numpy.interp(a, a_mod, res) # pyright:ignore[reportUnknownArgumentType] # re-sampled back to original timestamps return numpy.array(res).astype(numpy.double) def smooth_list(self, aa:'npt.NDArray[numpy.double]|npt.NDArray[numpy.floating]|Sequence[float]', b:'npt.NDArray[numpy.double]|npt.NDArray[numpy.floating]|Sequence[float]', window_len:int = 7, window:str = 'hanning', decay_weights:list[int]|None = None, decay_smoothing:bool = False, fromIndex:int = -1, toIndex:int = 0, re_sample:bool = True, back_sample:bool = True, a_lin:'npt.NDArray[numpy.double]|None' = None, delta:bool=False) -> 'npt.NDArray[numpy.double]': if len(aa) > 1 and len(aa) == len(b) and (self.filterDropOuts or window_len>2): #pylint: disable=E1103 # 1. truncate if fromIndex > -1: # if fromIndex is set, replace prefix up to fromIndex by None if toIndex==0: # no limit toIndex=len(aa) else: # smooth list on full length fromIndex = 0 toIndex = len(aa) a = numpy.array(aa[fromIndex:toIndex], dtype=numpy.double) # we mask the error value -1 and Numpy in the temperature array mb:numpy.ndarray[tuple[Literal[1]],numpy.dtype[numpy.float64]] = cast(numpy.ndarray[tuple[Literal[1]],numpy.dtype[numpy.float64]], numpy.ma.masked_equal(b[fromIndex:toIndex], -1)) # split in masked and unmasked_slices = [(x,False) for x in numpy.ma.clump_unmasked(mb)] # type:ignore[no-untyped-call] # ty:ignore[ignore] # the valid readings masked_slices = [(x,True) for x in numpy.ma.clump_masked(mb)] # type:ignore[no-untyped-call] # ty:ignore[ignore] # the dropped values sorted_slices = sorted(unmasked_slices + masked_slices, key=lambda tup: tup[0].start) # pyright:ignore[reportUnknownArgumentType] # pyright: ignore[reportGeneralTypeIssues] b_smoothed:list[npt.NDArray[numpy.double]] = [] # pyright:ignore[reportUnknownArgumentType] # b_smoothed collects the smoothed segments in order b_smoothed.append(numpy.full(fromIndex, numpy.nan, dtype=numpy.double)) # pyright:ignore[reportUnknownArgumentType] # append initial segment to the list of resulting segments # we just smooth the unmsked slices and add the unmasked slices with NaN values for (s, m) in sorted_slices: if m: # a slice with all masked (invalid) readings b_smoothed.append(numpy.full(s.stop - s.start, numpy.nan, dtype=numpy.double)) # pyright:ignore[reportUnknownArgumentType] else: # a slice with proper data b_smoothed.append(self.smooth_slice(a[s], mb[s], window_len, window, decay_weights, decay_smoothing, re_sample, back_sample, a_lin, delta)) # pyright:ignore[reportUnknownArgumentType] b_smoothed.append(numpy.full(len(a)-toIndex, numpy.nan, dtype=numpy.double)) # append the final segment to the list of resulting segments return numpy.concatenate(b_smoothed) bb = numpy.array(b, dtype=numpy.double) bb[bb == -1] = numpy.nan return bb @staticmethod def medfilt(x:'npt.NDArray[numpy.double]', k:int) -> 'npt.NDArray[numpy.double]': """Apply a length-k median filter to a 1D array x. Boundaries are extended by repeating endpoints. """ assert k % 2 == 1, 'Median filter length must be odd.' assert x.ndim == 1, 'Input must be one-dimensional.' if len(x) == 0: return x k2 = (k - 1) // 2 y = numpy.zeros ((len (x), k), dtype=x.dtype) y[:,k2] = x for i in range (k2): j = k2 - i y[j:,i] = x[:-j] y[:j,i] = x[0] y[:-j,-(i+1)] = x[j:] y[-j:,-(i+1)] = x[-1] return numpy.median(y, axis=1) @staticmethod def arrayRoR(tx:'npt.NDArray[numpy.double]', temp:'npt.NDArray[numpy.double]', wsize:int) -> 'npt.NDArray[numpy.floating]': # with wsize >=1 # length compensation done downstream, not necessary here! with warnings.catch_warnings(): # suppress warning if time difference is 0 which leads to a div by zero resulting in a warning and an inf value warnings.simplefilter('ignore') return (temp[wsize:] - temp[:-wsize]) / ((tx[wsize:] - tx[:-wsize])/60.) @staticmethod def polyRoR(tx:'npt.NDArray[numpy.double]', temp:'npt.NDArray[numpy.double]', wsize:int, i:int) -> float: if i == 0: # we duplicate the first possible RoR value instead of returning a 0 i = 1 if 0 < i < min(len(tx), len(temp)): with warnings.catch_warnings(): warnings.simplefilter('ignore') left_index = max(0,i-wsize) LS_fit = numpy.polynomial.polynomial.polyfit(tx[left_index:i+1],temp[left_index:i+1], 1) return float(LS_fit[1]*60.) else: return 0 def updateProjection(self) -> None: try: # projections are only drawn after CHARGE and before DROP if self.ax is not None and self.timeindex[0] != -1 and self.timeindex[6] == 0 and len(self.timex)>0: charge = self.timex[self.timeindex[0]] # in data time (corresponds to display time 00:00) now = self.timex[-1] # in data time (incl. time to charge) _,xlim_right = self.ax.get_xlim() # in data time like timex (incl. time to charge) if self.projectionmode == 0 or (self.projectionmode == 1 and (self.timex[-1]-charge)<=60*5): # linear temperature projection mode based on current RoR #calculate the temperature endpoint at endofx according to the latest rate of change if self.l_BTprojection is not None: if self.BTcurve and len(self.unfiltereddelta2_pure) > 0 and len(self.ctemp2) > 0 and self.ctemp2[-1] is not None and self.ctemp2[-1] != -1 and not math.isnan(self.ctemp2[-1]): # projection extended to the plots current endofx left = now right = max(left, xlim_right + charge) # never have the right point be left of left;) BTprojection = self.ctemp2[-1] + self.unfiltereddelta2_pure[-1]*(right - left)/60. # pyrefly: ignore[unsupported-operation] #plot projection self.BTprojection_tx = [left,right] self.BTprojection_temp = [self.ctemp2[-1], BTprojection] # pyrefly: ignore[bad-assignment] else: self.BTprojection_tx = [] self.BTprojection_temp = [] self.l_BTprojection.set_data(self.BTprojection_tx, self.BTprojection_temp) if self.l_ETprojection is not None: if self.ETcurve and len(self.unfiltereddelta1_pure) > 0 and len(self.ctemp1) > 0 and self.ctemp1[-1] is not None and self.ctemp1[-1] != -1 and not math.isnan(self.ctemp1[-1]): # projection extended to the plots current endofx left = now right = max(left,xlim_right + charge) # never have the right point be left of left;) ETprojection = self.ctemp1[-1] + self.unfiltereddelta1_pure[-1]*(right - left)/60. # pyrefly: ignore[unsupported-operation] #plot projection self.ETprojection_tx = [left,right] self.ETprojection_temp = [self.ctemp1[-1], ETprojection] # pyrefly: ignore[bad-assignment] else: self.ETprojection_tx = [] self.ETprojection_temp = [] self.l_ETprojection.set_data(self.ETprojection_tx, self.ETprojection_temp) # quadratic temperature projection based on linear RoR approximation # only active 5min after CHARGE elif self.projectionmode == 1 and (self.timex[-1]-charge)>60*5: delta_interval_BT = max(10, self.deltaBTsamples) # at least a span of 10 readings delta_interval_ET = max(10, self.deltaETsamples) # at least a span of 10 readings deltadeltalimit = 0.002 delay = self.delay/1000. # NOTE: we use the unfiltered deltas here to make this work also with a delta symbolic formula like x/2 to render RoR in C/30sec if self.l_BTprojection is not None: if (len(self.ctemp2) > 0 and self.ctemp2[-1] is not None and self.ctemp2[-1] != -1 and not math.isnan(self.ctemp2[-1]) and len(self.unfiltereddelta2_pure)>delta_interval_BT and self.unfiltereddelta2_pure[-1] and self.unfiltereddelta2_pure[-1]>0 and self.unfiltereddelta2_pure[-delta_interval_BT] and self.unfiltereddelta2_pure[-delta_interval_BT]>0): deltadelta_secsec = (((self.unfiltereddelta2_pure[-1] - self.unfiltereddelta2_pure[-delta_interval_BT])/60) / (now - self.timex[-delta_interval_BT])) # linear BT RoRoR C/sec/sec # limit deltadelta deltadelta_secsec = max(-deltadeltalimit,min(deltadeltalimit,deltadelta_secsec)) xpoints = numpy.arange(now, xlim_right, delay) ypoints = [self.ctemp2[-1]] delta_sec = self.unfiltereddelta2_pure[-1]/60 for _ in range(len(xpoints)-1): ypoints.append(ypoints[-1] + delta_sec*delay) # pyrefly: ignore[unsupported-operation] delta_sec = delta_sec + deltadelta_secsec*delay #plot BT curve self.BTprojection_tx = xpoints.tolist() self.BTprojection_temp = ypoints # pyrefly: ignore[bad-assignment] else: self.BTprojection_tx = [] self.BTprojection_temp = [] self.l_BTprojection.set_data(self.BTprojection_tx, self.BTprojection_temp) if self.l_ETprojection is not None: if (len(self.ctemp1) > 0 and self.ctemp1[-1] is not None and self.ctemp1[-1] != -1 and not math.isnan(self.ctemp1[-1]) and len(self.unfiltereddelta1_pure)>delta_interval_BT and self.unfiltereddelta1_pure[-1] and self.unfiltereddelta1_pure[-1]>0 and self.unfiltereddelta1_pure[-delta_interval_BT] and self.unfiltereddelta1_pure[-delta_interval_BT]>0): deltadelta_secsec = (((self.unfiltereddelta1_pure[-1] - self.unfiltereddelta1_pure[-delta_interval_ET])/60) / (self.timex[-1] - self.timex[-delta_interval_ET])) # linear ET RoRoR C/sec/sec # limit deltadelta deltadelta_secsec = max(-deltadeltalimit,min(deltadeltalimit,deltadelta_secsec)) xpoints = numpy.arange(now, xlim_right, delay) ypoints = [self.ctemp1[-1]] delta_sec = self.unfiltereddelta1_pure[-1]/60 for _ in range(len(xpoints)-1): ypoints.append(ypoints[-1] + delta_sec*delay) # pyrefly: ignore[unsupported-operation] delta_sec = delta_sec + deltadelta_secsec*delay #plot ET curve self.ETprojection_tx = xpoints.tolist() self.ETprojection_temp = ypoints # pyrefly: ignore[bad-assignment] else: self.ETprojection_tx = [] self.ETprojection_temp = [] self.l_ETprojection.set_data(self.ETprojection_tx, self.ETprojection_temp) # RoR projections if self.projectDeltaFlag and (self.timex[-1]-charge)>60*5: delay = self.delay/1000. if self.l_DeltaBTprojection is not None: delta_interval_BT = max(10, self.deltaBTsamples*2) # at least a span of 10 readings if (self.DeltaBTflag and len(self.delta2)>0 and len(self.delta2)>delta_interval_BT): d2_last = self.delta2[-1] d2_left = self.delta2[-delta_interval_BT] if d2_last is not None and d2_left is not None and d2_last>0 and d2_left>0: # compute deltadelta_secsec from delta2 adjusted to delta math formulas deltadelta_secsec = (((d2_last - d2_left)/60) / (now - self.timex[-delta_interval_BT])) # linear BT RoRoR C/sec/sec left = now right = max(left, xlim_right) # never have the right point be left of left;) DeltaBTprojection = d2_last + deltadelta_secsec * (right - left) * 60 # projection extended to the plots current endofx self.DeltaBTprojection_tx = [left,right] self.DeltaBTprojection_temp = [d2_last, DeltaBTprojection] else: self.DeltaBTprojection_tx = [] self.DeltaBTprojection_temp = [] else: self.DeltaBTprojection_tx = [] self.DeltaBTprojection_temp = [] self.l_DeltaBTprojection.set_data(self.DeltaBTprojection_tx, self.DeltaBTprojection_temp) if self.l_DeltaETprojection is not None: delta_interval_ET = max(10, self.deltaETsamples*2) # at least a span of 10 readings if (self.DeltaETflag and len(self.delta1)>0 and len(self.delta1)>delta_interval_ET): d1_last = self.delta1[-1] d1_left = self.delta1[-delta_interval_ET] if d1_last is not None and d1_left is not None and d1_last>0 and d1_left>0: # compute deltadelta_secsec from delta1 adjusted to delta math formulas deltadelta_secsec = (((d1_last - d1_left)/60) / (now - self.timex[-delta_interval_ET])) # linear ET RoRoR C/sec/sec left = now right = max(left, xlim_right) # never have the right point be left of left;) DeltaETprojection = d1_last + deltadelta_secsec * (right - left) * 60 # projection extended to the plots current endofx self.DeltaETprojection_tx = [left,right] self.DeltaETprojection_temp = [d1_last, DeltaETprojection] else: self.DeltaETprojection_tx = [] self.DeltaETprojection_temp = [] else: self.DeltaETprojection_tx = [] self.DeltaETprojection_temp = [] self.l_DeltaETprojection.set_data(self.DeltaETprojection_tx, self.DeltaETprojection_temp) else: self.BTprojection_tx = [] self.BTprojection_temp = [] self.ETprojection_tx = [] self.ETprojection_temp = [] self.DeltaBTprojection_tx = [] self.DeltaBTprojection_temp = [] self.DeltaETprojection_tx = [] self.DeltaETprojection_temp = [] if self.l_BTprojection is not None: self.l_BTprojection.set_data([],[]) if self.l_ETprojection is not None: self.l_ETprojection.set_data([],[]) if self.l_DeltaBTprojection is not None: self.l_DeltaBTprojection.set_data([],[]) if self.l_DeltaETprojection is not None: self.l_DeltaETprojection.set_data([],[]) except Exception as ex: # pylint: disable=broad-except _log.exception(ex) _, _, exc_tb = sys.exc_info() self.adderror((QApplication.translate('Error Message','Exception:') + ' updateProjection() {0}').format(str(ex)),getattr(exc_tb, 'tb_lineno', '?')) if self.l_BTprojection is not None: self.l_BTprojection.set_data([],[]) if self.l_ETprojection is not None: self.l_ETprojection.set_data([],[]) if self.l_DeltaBTprojection is not None: self.l_DeltaBTprojection.set_data([],[]) if self.l_DeltaETprojection is not None: self.l_DeltaETprojection.set_data([],[]) def process_values(self,tx,t1,t2): local_flagstart=self.flagstart sample_timex=self.timex sample_temp1=self.temp1;sample_temp2=self.temp2 sample_ctimex1=self.ctimex1;sample_ctimex2=self.ctimex2 sample_ctemp1=self.ctemp1;sample_ctemp2=self.ctemp2 sample_tstemp1=self.tstemp1;sample_tstemp2=self.tstemp2 sample_unfiltereddelta1=self.unfiltereddelta1;sample_unfiltereddelta2=self.unfiltereddelta2 sample_delta1=self.delta1;sample_delta2=self.delta2 t1_prev = t2_prev = None t1_prevprev = t2_prevprev = None if (self.minmaxLimits or self.dropSpikes or self.dropDuplicates): if len(sample_temp1)>0: t1_prev = sample_temp1[-1] if len(sample_temp1)>1: t1_prevprev = sample_temp1[-2] if len(sample_temp2)>0: t2_prev = sample_temp2[-1] if len(sample_temp2)>1: t2_prevprev = sample_temp2[-2] if not self.dummy_or_special_device(self.device, 0): t1 = self.inputFilter(sample_timex,sample_temp1,tx,t1) if not self.dummy_or_special_device(self.device, 0): t2 = self.inputFilter(sample_timex,sample_temp2,tx,t2,True) length_of_qmc_timex = len(sample_timex) # now copy the destructively modified values from temp1/2 to ctemp1/2 if any (to ensure to pick the right elements we compare the timestamps at those indices) if (self.minmaxLimits or self.dropSpikes or self.dropDuplicates): if len(sample_ctimex1)>0: if t1_prev is not None and sample_ctimex1[-1] == sample_timex[-1] and t1_prev != sample_temp1[-1]: sample_ctemp1[-1] = sample_temp1[-1] if len(sample_ctimex1)>1 and t1_prevprev is not None and sample_ctimex1[-2] == sample_timex[-2] and t1_prevprev != sample_temp1[-2]: sample_ctemp1[-2] = sample_temp1[-2] if len(sample_ctimex2)>0: if t2_prev is not None and sample_ctimex2[-1] == sample_timex[-1] and t2_prev != sample_temp2[-1]: sample_ctemp2[-1] = sample_temp2[-1] if len(sample_ctimex2)>1 and t2_prevprev is not None and sample_ctimex2[-2] == sample_timex[-2] and t2_prevprev != sample_temp2[-2]: sample_ctemp2[-2] = sample_temp2[-2] t1_final = t1 t2_final = t2 sample_temp1.append(t1_final) sample_temp2.append(t2_final) sample_timex.append(tx) length_of_qmc_timex += 1 if t1_final != -1: sample_ctimex1.append(tx) sample_ctemp1.append(t1_final) elif len(sample_temp1)>(self.interpolatemax+1) and all(v == -1 for v in sample_temp1[-(self.interpolatemax+1):]): sample_ctimex1.append(tx) sample_ctemp1.append(None) if t2_final != -1: sample_ctimex2.append(tx) sample_ctemp2.append(t2_final) elif len(sample_temp2)>(self.interpolatemax+1) and all(v == -1 for v in sample_temp2[-(self.interpolatemax+1):]): sample_ctimex2.append(tx) sample_ctemp2.append(None) #we populate the temporary smoothed ET/BT data arrays (with readings cleansed from -1 dropouts) cf = self.curvefilter if self.temp_decay_weights is None or len(self.temp_decay_weights) != cf: # recompute only on changes self.temp_decay_weights = list(numpy.arange(1,cf+1)) # we don't smooth st'x if last, or butlast temperature value were a drop-out not to confuse the RoR calculation if -1 in sample_temp1[-(cf+1):]: dw1 = [1] else: dw1 = self.temp_decay_weights if -1 in sample_temp2[-(cf+1):]: dw2 = [1] else: dw2 = self.temp_decay_weights # average smoothing if len(sample_ctemp1) > 0: st1 = self.decay_average(sample_ctimex1,sample_ctemp1,dw1) else: st1 = -1 if len(sample_ctemp2) > 0: st2 = self.decay_average(sample_ctimex2,sample_ctemp2,dw2) else: st2 = -1 # # we apply a minimal live median spike filter # if self.filterDropOuts: # if st1 is not None and st1 != -1: # st1 = self.liveMedianETfilter(st1) # if st1 is not None and st2 != -1: # st2 = self.liveMedianBTfilter(st2) # register smoothed values sample_tstemp1.append(st1) sample_tstemp2.append(st2) rateofchange1plot:float|None rateofchange2plot:float|None #we need a minimum of two readings to calculate rate of change if length_of_qmc_timex > 1: # compute T1 RoR try: if t1_final == -1 or len(sample_ctimex1)<2: # we repeat the last RoR if underlying temperature dropped if sample_unfiltereddelta1: self.rateofchange1 = sample_unfiltereddelta1[-1] else: self.rateofchange1 = 0. else: # normal data received # Delta T = (changeTemp/ChangeTime)*60. = degrees per minute; left_index = min(len(sample_ctimex1),len(sample_tstemp1),max(2, self.deltaETsamples + 1)) # ****** Instead of basing the estimate on the window extremal points, # grab the full set of points and do a formal LS solution to a straight line and use the slope estimate for RoR if self.polyfitRoRcalc: try: time_vec = sample_ctimex1[-left_index:] temp_samples = sample_tstemp1[-left_index:] with warnings.catch_warnings(): warnings.simplefilter('ignore') # using stable polyfit from numpy polyfit module LS_fit = numpy.polynomial.polynomial.polyfit(time_vec, temp_samples, 1) self.rateofchange1 = LS_fit[1]*60. except Exception: # pylint: disable=broad-except # a numpy/OpenBLAS polyfit bug can cause polyfit to throw an exception "SVD did not converge in Linear Least Squares" on Windows Windows 10 update 2004 # https://github.com/numpy/numpy/issues/16744 # we fall back to the two point algo timed = sample_ctimex1[-1] - sample_ctimex1[-left_index] #time difference between last self.deltaETsamples readings self.rateofchange1 = ((sample_tstemp1[-1] - sample_tstemp1[-left_index])/timed)*60. #delta ET (degrees/minute) else: timed = sample_ctimex1[-1] - sample_ctimex1[-left_index] #time difference between last self.deltaETsamples readings self.rateofchange1 = ((sample_tstemp1[-1] - sample_tstemp1[-left_index])/timed)*60. #delta ET (degrees/minute) except Exception as e: # pylint: disable=broad-except _log.error(e) self.rateofchange1 = 0. # compute T2 RoR try: if t2_final == -1 or len(sample_ctimex2)<2: # we repeat the last RoR if underlying temperature dropped if sample_unfiltereddelta2: self.rateofchange2 = sample_unfiltereddelta2[-1] else: self.rateofchange2 = 0. else: # normal data received # Delta T = (changeTemp/ChangeTime)*60. = degrees per minute; left_index = min(len(sample_ctimex2),len(sample_tstemp2),max(2, self.deltaBTsamples + 1)) # ****** Instead of basing the estimate on the window extremal points, # grab the full set of points and do a formal LS solution to a straight line and use the slope estimate for RoR if self.polyfitRoRcalc: try: time_vec = sample_ctimex2[-left_index:] temp_samples = sample_tstemp2[-left_index:] with warnings.catch_warnings(): warnings.simplefilter('ignore') LS_fit = numpy.polynomial.polynomial.polyfit(time_vec, temp_samples, 1) self.rateofchange2 = LS_fit[1]*60. except Exception: # pylint: disable=broad-except # a numpy/OpenBLAS polyfit bug can cause polyfit to throw an exception "SVD did not converge in Linear Least Squares" on Windows Windows 10 update 2004 # https://github.com/numpy/numpy/issues/16744 # we fall back to the two point algo timed = sample_ctimex2[-1] - sample_ctimex2[-left_index] #time difference between last self.deltaBTsamples readings self.rateofchange2 = ((sample_tstemp2[-1] - sample_tstemp2[-left_index])/timed)*60. #delta BT (degrees/minute) else: timed = sample_ctimex2[-1] - sample_ctimex2[-left_index] #time difference between last self.deltaBTsamples readings self.rateofchange2 = ((sample_tstemp2[-1] - sample_tstemp2[-left_index])/timed)*60. #delta BT (degrees/minute) except Exception as e: # pylint: disable=broad-except _log.error(e) self.rateofchange1 = 0. # self.unfiltereddelta{1,2}_pure contain the RoR values respecting the delta_span, but without any delta smoothing NOR delta mathformulas applied self.unfiltereddelta1_pure.append(self.rateofchange1) self.unfiltereddelta2_pure.append(self.rateofchange2) # apply the math formula before the delta smoothing if self.DeltaETfunction.strip(): try: self.rateofchange1 = self.eval_math_expression(self.DeltaETfunction,tx,RTsname='R1',RTsval=self.rateofchange1) except Exception as e: # pylint: disable=broad-except _log.exception(e) if self.DeltaBTfunction.strip(): try: self.rateofchange2 = self.eval_math_expression(self.DeltaBTfunction,tx,RTsname='R2',RTsval=self.rateofchange2) except Exception as e: # pylint: disable=broad-except _log.exception(e) # RoR Smooth Spikes (Drop Filter, Median Filter) is only applied offline # # unfiltereddelta1/2 contains the RoRs respecting the delta_span, but without any delta smoothing AND delta mathformulas applied # # we apply a minimal live median spike filter minimizing the delay by choosing a window smaller than in the offline medfilt # if self.filterDropOuts: # if self.rateofchange1 is not None and self.rateofchange1 != -1: # self.rateofchange1 = self.liveMedianETRoRfilter(self.rateofchange1) # if self.rateofchange2 is not None and self.rateofchange2 != -1: # self.rateofchange2 = self.liveMedianBTRoRfilter(self.rateofchange2) sample_unfiltereddelta1.append(self.rateofchange1) sample_unfiltereddelta2.append(self.rateofchange2) ####### filter deltaBT deltaET # decay smoothing if self.deltaETfilter: user_filter = int(round(self.deltaETfilter/2.)) if user_filter and length_of_qmc_timex > user_filter and (len(sample_unfiltereddelta1) > user_filter): if self.decay_weights is None or len(self.decay_weights) != user_filter: # recompute only on changes self.decay_weights = list(numpy.arange(1,user_filter+1)) self.rateofchange1 = self.decay_average(sample_timex,sample_unfiltereddelta1,self.decay_weights) if self.deltaBTfilter: user_filter = int(round(self.deltaBTfilter/2.)) if user_filter and length_of_qmc_timex > user_filter and (len(sample_unfiltereddelta2) > user_filter): if self.decay_weights is None or len(self.decay_weights) != user_filter: # recompute only on changes self.decay_weights = list(numpy.arange(1,user_filter+1)) self.rateofchange2 = self.decay_average(sample_timex,sample_unfiltereddelta2,self.decay_weights) rateofchange1plot = self.rateofchange1 rateofchange2plot = self.rateofchange2 else: sample_unfiltereddelta1.append(0.) sample_unfiltereddelta2.append(0.) self.rateofchange1,self.rateofchange2,rateofchange1plot,rateofchange2plot = 0.,0.,0.,0. # limit displayed RoR #(only before TP is recognized) # WHY? if self.RoRlimitFlag: # not self.TPalarmtimeindex and self.RoRlimitFlag: if not max(-self.maxRoRlimit,self.RoRlimitm) < rateofchange1plot < min(self.maxRoRlimit,self.RoRlimit): rateofchange1plot = None if not max(-self.maxRoRlimit,self.RoRlimitm) < rateofchange2plot < min(self.maxRoRlimit,self.RoRlimit): rateofchange2plot = None # append new data to the rateofchange arrays sample_delta1.append(rateofchange1plot) sample_delta2.append(rateofchange2plot)