Compute a baseline trend of a signal
430082
def detrend(arr, x=None, deg=5, tol=1e-3, maxloop=10):
    """"""

    xx = numpy.arange(len(arr)) if x is None else x
    base = arr.copy()
    trend = base
    pol = numpy.ones((deg + 1,))
    for _ in range(maxloop):
        pol_new = numpy.polyfit(xx, base, deg)
        pol_norm = numpy.linalg.norm(pol)
        diff_pol_norm = numpy.linalg.norm(pol - pol_new)
        if diff_pol_norm / pol_norm < tol:
            break
        pol = pol_new
        trend = numpy.polyval(pol, xx)
        base = numpy.minimum(base, trend)
    return trend
936468
def flag_calcmad(data):
    

    absdev = np.abs(data - np.median(data))
    return np.median(absdev)
934657
def _rm_outlier_by_amp(self, params, model, signal, ii):
        
        maxamps = np.nanmax(np.abs(model),0)
        z_score = (maxamps - np.nanmean(maxamps,0))/np.nanstd(maxamps,0)
        with warnings.catch_warnings():
            warnings.simplefilter("ignore")
            outlier_idx = np.where(np.abs(z_score)>2.0)[0]
            nan_idx = np.where(np.isnan(params))[0]
            outlier_idx = np.unique(np.hstack([nan_idx, outlier_idx]))
            ii[outlier_idx] = 0
            model[outlier_idx] = np.nan
            signal[outlier_idx] = np.nan
            params[outlier_idx] = np.nan

        return model, signal, params, ii
671999
def _update_noise(self, peak_num):
        
        i = self.peak_inds_i[peak_num]
        self.noise_amp_recent = (0.875*self.noise_amp_recent
                                 + 0.125*self.sig_i[i])
        return
347171
def _get_noise(self, residuals: np.ndarray) -> float:
        
        return np.mean(np.abs(residuals))
571831
def median(self):
        """"""
        mu = self.mean()
        ret_val = math.exp(mu)
        if math.isnan(ret_val):
            ret_val = float("inf")
        return ret_val
189850
def estimate_cp_pval(ts, method="mean"):
    

    """ 
        ts: time series
        method: look for a single changepoint in 'mean' , 'var', 'mean and var'
        Returns: returns index of the changepoint, and pvalue. Here pvalue = 1
                 means statistically significant
    """
    robjects.r("library(changepoint)")
    method_map = {
        "mean": "cpt.mean({}, class=FALSE)",
        "var": "cpt.var({}, class=FALSE)",
        "meanvar": "cpt.meanvar({}, class=FALSE)",
    }
    mt = robjects.FloatVector(ts)
    robjects.globalenv["mt"] = mt
    cmd = method_map[method].format("mt")
    robjects.globalenv["mycpt"] = robjects.r(cmd)
    ecp_pval = robjects.r("mycpt")
    ecp = ecp_pval[0]
    pval = ecp_pval[1]
    return ecp, pval
725164
def median(series):
    

    if np.issubdtype(series.dtype, np.number):
        return series.median()
    else:
        return np.nan
1089715
def ystep(self):
        
        """

        rho = self.Nb * self.rho
        mAXU = np.mean(self.AX + self.U, axis=-1)
        self.Y[:] = self.prox_g(mAXU, rho)
337141
def Kdiag(self, X, target):
        """"""
        ly=1/self.lengthscaleY
        lu=np.sqrt(3)/self.lengthscaleU
        #ly=self.lengthscaleY
        #lu=self.lengthscaleU
        
        k1 = (2*lu+ly)/(lu+ly)**2
        k2 = (ly-2*lu + 2*lu-ly ) / (ly-lu)**2 
        k3 = 1/(lu+ly) + (lu)/(lu+ly)**2 

        np.add(self.varianceU*self.varianceY*(k1+k2+k3), target, target)
254073
def moving_average_convergence(data, nslow=26, nfast=12):
    
    emaslow = moving_average(data, nslow, type='exponential')
    emafast = moving_average(data, nfast, type='exponential')
    return emaslow, emafast, emafast - emaslow
1372395
def loadCurve(data, groups, thresholds, absvals, fs, xlabels):
        """"""
        xlims = (xlabels[0], xlabels[-1])
        pw = ProgressWidget(groups, xlims)
        spike_counts = []
        # skip control
        for itrace in range(data.shape[0]):
            count = 0
            for ichan in range(data.shape[2]):
                flat_reps = data[itrace,:,ichan,:].flatten()
                count += len(spikestats.spike_times(flat_reps, thresholds[ichan], fs, absvals[ichan]))
            spike_counts.append(count/(data.shape[1]*data.shape[2])) #mean spikes per rep

        i = 0
        for g in groups:
            for x in xlabels:
                pw.setPoint(x, g, spike_counts[i])
                i +=1

        return pw
174742
def simple_moving_average(data, period):
    
    catch_errors.check_for_period_error(data, period)
    # Mean of Empty Slice RuntimeWarning doesn't affect output so it is
    # supressed
    with warnings.catch_warnings():
        warnings.simplefilter("ignore", category=RuntimeWarning)
        sma = [np.mean(data[idx-(period-1):idx+1]) for idx in range(0, len(data))]
    sma = fill_for_noncomputable_vals(data, sma)
    return sma
444484
def Kdiag(self, X):
        """"""
        K_ff = np.zeros(X.shape[0])
        for i,x in enumerate(X):
            K_ff[i] = self.k_ff(x[0],x[0],self.lengthscale[0])
        return K_ff * self.variances[0]
226413
def plot_signal(signal, sig_len, n_sig, fs, time_units, sig_style, axes):
    ""

    # Extend signal style if necesary
    if len(sig_style) == 1:
        sig_style = n_sig * sig_style

    # Figure out time indices
    if time_units == 'samples':
        t = np.linspace(0, sig_len-1, sig_len)
    else:
        downsample_factor = {'seconds':fs, 'minutes':fs * 60,
                             'hours':fs * 3600}
        t = np.linspace(0, sig_len-1, sig_len) / downsample_factor[time_units]

    # Plot the signals
    if signal.ndim == 1:
        axes[0].plot(t, signal, sig_style[0], zorder=3)
    else:
        for ch in range(n_sig):
            axes[ch].plot(t, signal[:,ch], sig_style[ch], zorder=3)
718579
def llinregress(x,y):
    
    TINY = 1.0e-20
    if len(x) != len(y):
        raise ValueError('Input values not paired in linregress.  Aborting.')
    n = len(x)
    x = [float(_) for _ in x]
    y = [float(_) for _ in y]
    xmean = mean(x)
    ymean = mean(y)
    r_num = float(n*(summult(x,y)) - sum(x)*sum(y))
    r_den = math.sqrt((n*ss(x) - square_of_sums(x))*(n*ss(y)-square_of_sums(y)))
    r = r_num / r_den
    z = 0.5*math.log((1.0+r+TINY)/(1.0-r+TINY))
    df = n-2
    t = r*math.sqrt(df/((1.0-r+TINY)*(1.0+r+TINY)))
    prob = betai(0.5*df,0.5,df/(df+t*t))
    slope = r_num / float(n*ss(x) - square_of_sums(x))
    intercept = ymean - slope*xmean
    sterrest = math.sqrt(1-r*r)*samplestdev(y)
    return slope, intercept, r, prob, sterrest
847218
function(signal1, signal2, options){
    options = options || {};
    options.epsilon = options.epsilon || 1E-08;
    var sig1 = Signal(signal1);
    var sig2 = Signal(signal2);
    if(signal1.length != signal2.length) return false;
    var diffSQ = _.reduce(_.range(signal1.length), function(d,idx){
      var diff = (sig1[idx].value-sig2[idx].value);
        return d + diff * diff;
    },0);
    return diffSQ/signal1.length < options.epsilon;
  }