Airline delay data used in Gaussian Processes for Big Data by Hensman, Fusi and Lawrence
190528
def airline_delay(data_set='airline_delay', num_train=700000, num_test=100000, seed=default_seed):
    """"""

    if not data_available(data_set):
        download_data(data_set)

    dir_path = os.path.join(data_path, data_set)
    filename = os.path.join(dir_path, 'filtered_data.pickle')

    # 1. Load the dataset
    import pandas as pd
    data = pd.read_pickle(filename)

    # WARNING: removing year
    data.pop('Year')

    # Get data matrices
    Yall = data.pop('ArrDelay').values[:,None]
    Xall = data.values

    # Subset the data (memory!!)
    all_data = num_train+num_test
    Xall = Xall[:all_data]
    Yall = Yall[:all_data]

    # Get testing points
    np.random.seed(seed=seed)
    N_shuffled = permute(Yall.shape[0])
    train, test = N_shuffled[num_test:], N_shuffled[:num_test]
    X, Y = Xall[train], Yall[train]
    Xtest, Ytest = Xall[test], Yall[test]
    covariates =  ['month', 'day of month', 'day of week', 'departure time', 'arrival time', 'air time', 'distance to travel', 'age of aircraft / years']
    response = ['delay']
    return data_details_return({'X': X, 'Y': Y, 'Xtest': Xtest, 'Ytest': Ytest, 'seed' : seed, 'info': "Airline delay data used for demonstrating Gaussian processes for big data.", 'covariates': covariates, 'response': response}, data_set)
37125
def _hotellings_time(self):
        
        # get residuals
        df = self._df
        R = self._A - np.dot(self._X, self._Bhat)
        Sigma_Z = np.dot(R.T,R)*(1./df)
        # compute p-values
        T_2_list = []
        p_value_list = []
        for i in np.arange(0,self._Bhat.shape[0]):
            Bstar       = self._Bhat[i,np.arange(0,np.shape(self._A)[1])]
            lstar       = float(np.shape(Bstar)[0])
            cx          = self._Cx[i,i]
            Einv        = np.linalg.pinv(Sigma_Z)
            Zs          = Bstar/np.sqrt(cx)
            T_2         = ((df - lstar + 1.)/(df*lstar)*np.dot(np.dot(Zs,Einv),Zs.T))
            p_value     = 1. - stats.f.cdf(T_2, lstar, df - lstar + 1.)

            p_value_list.append(p_value)
            T_2_list.append(T_2)
        # save pvalue results in a table
        self._results = [['Comparison','Hotellings T^2', "p-value"]]
        for i in np.arange(0,len(self._col_names)):
            self._results.append([self._col_names[i], T_2_list[i], p_value_list[i]])
523794
def map_predict(interface, state, label, inp):
    
    import numpy as np
    out = interface.output(0)

    continuous = [j for i, j in enumerate(state["X_indices"]) if
                  state["X_meta"][i] == "c"]  # indices of continuous features
    discrete = [j for i, j in enumerate(state["X_indices"]) if
                state["X_meta"][i] == "d"]  # indices of discrete features

    cont = True if len(continuous) > 0 else False  # enables calculation of Gaussian probabilities
    disc = True if len(discrete) > 0 else False  # enables calculation of multinomial probabilities.

    for row in inp:
        row = row.strip().split(state["delimiter"])
        if len(row) > 1:  # if row is empty
            # set id of a sample
            x_id = "" if state["id_index"] == -1 else row[state["id_index"]]
            # initialize prior probability for all labels
            probs = state["fit_model"]["prior_log"]

            if cont:  # continuous features
                x = np.array([(0 if row[j] in state["missing_vals"] else float(row[j])) for j in
                              continuous])  # sets selected features of the sample
                # Gaussian distribution
                probs = probs - 0.5 * np.sum(
                    np.true_divide((x - state["fit_model"]["mean"]) ** 2, state["fit_model"]["var"]) +
                    state["fit_model"]["var_log"], axis=1)

            if disc:  # discrete features
                # multinomial distribution
                probs = probs + np.sum(
                    [(0 if row[i] in state["missing_vals"] else state["fit_model"].get((str(i), row[i]), np.zeros(1)))
                     for i in discrete], axis=0)

            # normalize by P(x) = P(f_1, ..., f_n)
            log_prob_x = np.log(np.sum(np.exp(probs)))
            probs = np.exp(np.array(probs) - log_prob_x)
            # Predicted label is the one with highest probability
            y_predicted = max(zip(probs, state["fit_model"]["y_labels"]))[1]
            out.add(x_id, (y_predicted, probs.tolist()))
944544
def lnlike(x, star):
    """"""
    ll = lnprior(x)
    if np.isinf(ll):
        return ll, (np.nan, np.nan)
    per, t0, b = x
    model = TransitModel('b', per=per, t0=t0, b=b, rhos=10.)(star.time)
    like, d, vard = star.lnlike(model, full_output=True)
    ll += like
    return ll, (d,)
712857
def delay(fn, opts, task, *args, **kwargs):
    
    if opts["stddev"] <= 0:
        t = opts["t"]
    elif opts["pdf"] == "gauss":
        t = random.gauss(opts["t"], opts["stddev"])
    elif opts["pdf"] == "uniform":
        t = random.uniform(opts["t"], opts["stddev"])
    else:
        raise ValueError("unknown delay decorator pdf '{}'".format(opts["pdf"]))

    time.sleep(t)

    return fn(task, *args, **kwargs)
494846
private void computeAxisVariance(List<P> points) {
		int numPoints = points.size();

		for( int i = 0; i < N; i++ ) {
			mean[i] = 0;
			var[i] = 0;
		}

		// compute the mean
		for( int i = 0; i < numPoints; i++ ) {
			P p = points.get(i);

			for( int j = 0; j < N; j++ ) {
				mean[j] += distance.valueAt(p,j);
			}
		}

		for( int i = 0; i < N; i++ ) {
			mean[i] /= numPoints;
		}

		// compute the variance * N
		for( int i = 0; i < numPoints; i++ ) {
			P p = points.get(i);

			for( int j = 0; j < N; j++ ) {
				double d = mean[j] - distance.valueAt(p,j);
				var[j] += d*d;
			}
		}
	}
1360060
def _compute_state_estimate(self):
        
        self.x.fill(0)
        for f, mu in zip(self.filters, self.mu):
            self.x += f.x * mu

        self.P.fill(0)
        for f, mu in zip(self.filters, self.mu):
            y = f.x - self.x
            self.P += mu * (outer(y, y) + f.P)
237590
def analyze(self, time: int=None) -> typing.Tuple[float, float, int, int]:
        
        if not time:
            time = int(len(self.times) / 2)
        if self.domains[time] == 0:
            yhat, residuals, mean_residual, error = self._get_fit(time)
            yhat_p = self.ddiff(yhat)
            yhat_pp = self.ddiff(yhat_p)
            noise = self._get_noise(residuals)
            curvature = (1 / self.ratio) * (1 / len(yhat_pp)) * np.sqrt(si.simps(yhat_pp ** 2))
            rng = (self.ratio * (np.max(self.averagedata[:, time]) -
                                 np.min(self.averagedata[:, time])))
            dmn = self.ratio * len(self.averagedata[:, time])

            self.noises[time] = np.log10(noise)
            self.curves[time] = np.log10(curvature)
            self.ranges[time] = np.log10(rng)
            self.domains[time] = np.log10(dmn)
        return self.noises[time], self.curves[time], self.ranges[time], self.domains[time]
1041939
private static int[] seed(SparseDataset data, int k) {
        int n = data.size();

        int[] y = new int[n];
        SparseArray centroid = data.get(Math.randomInt(n)).x;

        double[] D = new double[n];
        for (int i = 0; i < n; i++) {
            D[i] = Double.MAX_VALUE;
        }

        // pick the next center
        for (int i = 1; i < k; i++) {
            for (int j = 0; j < n; j++) {
                double dist = Math.JensenShannonDivergence(data.get(j).x, centroid);
                if (dist < D[j]) {
                    D[j] = dist;
                    y[j] = i - 1;
                }
            }

            double cutoff = Math.random() * Math.sum(D);
            double cost = 0.0;
            int index = 0;
            for (; index < n; index++) {
                cost += D[index];
                if (cost >= cutoff) {
                    break;
                }
            }

            centroid = data.get(index).x;
        }

        for (int j = 0; j < n; j++) {
            // compute the distance between this sample and the current center
            double dist = Math.JensenShannonDivergence(data.get(j).x, centroid);
            if (dist < D[j]) {
                D[j] = dist;
                y[j] = k - 1;
            }
        }

        return y;
    }
206592
def factorized_gaussian_noise(in_features, out_features, device):
    
    in_noise = scaled_noise(in_features, device=device)
    out_noise = scaled_noise(out_features, device=device)

    return out_noise.ger(in_noise), out_noise
139995
public static void recomputeCovarianceMatrices(Relation<? extends NumberVector> relation, WritableDataStore<double[]> probClusterIGivenX, List<? extends EMClusterModel<?>> models, double prior) {
    final int k = models.size();
    boolean needsTwoPass = false;
    for(EMClusterModel<?> m : models) {
      m.beginEStep();
      needsTwoPass |= m.needsTwoPass();
    }
    // First pass, only for two-pass models.
    if(needsTwoPass) {
      for(DBIDIter iditer = relation.iterDBIDs(); iditer.valid(); iditer.advance()) {
        double[] clusterProbabilities = probClusterIGivenX.get(iditer);
        NumberVector instance = relation.get(iditer);
        for(int i = 0; i < clusterProbabilities.length; i++) {
          final double prob = clusterProbabilities[i];
          if(prob > 1e-10) {
            models.get(i).firstPassE(instance, prob);
          }
        }
      }
      for(EMClusterModel<?> m : models) {
        m.finalizeFirstPassE();
      }
    }
    double[] wsum = new double[k];
    for(DBIDIter iditer = relation.iterDBIDs(); iditer.valid(); iditer.advance()) {
      double[] clusterProbabilities = probClusterIGivenX.get(iditer);
      NumberVector instance = relation.get(iditer);
      for(int i = 0; i < clusterProbabilities.length; i++) {
        final double prob = clusterProbabilities[i];
        if(prob > 1e-10) {
          models.get(i).updateE(instance, prob);
        }
        wsum[i] += prob;
      }
    }
    for(int i = 0; i < models.size(); i++) {
      // MLE / MAP
      final double weight = prior <= 0. ? wsum[i] / relation.size() : (wsum[i] + prior - 1) / (relation.size() + prior * k - k);
      models.get(i).finalizeEStep(weight, prior);
    }
  }
803198
def funcNrlTcMotPred(idxPrc,
                     varPixX,
                     varPixY,
                     NrlMdlChunk,
                     varNumTP,
                     aryBoxCar,  # aryCond
                     path,
                     varNumNrlMdls,
                     varNumMtDrctn,
                     varPar,
                     queOut):
    

#    # if hd5 method is used: open file for reading
#    filename = 'aryBoxCar' + str(idxPrc) + '.hdf5'
#    hdf5_path = os.path.join(path, filename)
#    fileH = tables.openFile(hdf5_path, mode='r')

    # Output array with pRF model time courses at all modelled standard
    # deviations for current pixel position:
    aryOut = np.empty((len(NrlMdlChunk), varNumTP, varNumMtDrctn),
                      dtype='float32')

    # Prepare status indicator if this is the first of the parallel processes:
    if idxPrc == 1:

        # We create a status indicator for the time consuming pRF model finding
        # algorithm. Number of steps of the status indicator:
        varStsStpSze = 20

        # Number of pRF models to fit:
        varNumLoops = varNumNrlMdls/varPar

        # Vector with pRF values at which to give status feedback:
        vecStatus = np.linspace(0,
                                varNumLoops,
                                num=(varStsStpSze+1),
                                endpoint=True)
        vecStatus = np.ceil(vecStatus)
        vecStatus = vecStatus.astype(int)

        # Vector with corresponding percentage values at which to give status
        # feedback:
        vecStatusPrc = np.linspace(0,
                                   100,
                                   num=(varStsStpSze+1),
                                   endpoint=True)
        vecStatusPrc = np.ceil(vecStatusPrc)
        vecStatusPrc = vecStatusPrc.astype(int)

        # Counter for status indicator:
        varCntSts01 = 0
        varCntSts02 = 0

    # Loop through all Gauss parameters that are in this chunk
    for idx, NrlMdlTrpl in enumerate(NrlMdlChunk):

        # Status indicator (only used in the first of the parallel
        # processes):
        if idxPrc == 1:

            # Status indicator:
            if varCntSts02 == vecStatus[varCntSts01]:

                # Prepare status message:
                strStsMsg = ('---------Progress: ' +
                             str(vecStatusPrc[varCntSts01]) +
                             ' % --- ' +
                             str(vecStatus[varCntSts01]) +
                             ' loops out of ' +
                             str(varNumLoops))

                print(strStsMsg)

                # Only increment counter if the last value has not been
                # reached yet:
                if varCntSts01 < varStsStpSze:
                    varCntSts01 = varCntSts01 + int(1)

        # x pos of Gauss model: NrlMdlTrpl[0]
        # y pos of Gauss model: NrlMdlTrpl[1]
        # std of Gauss model: NrlMdlTrpl[2]
        # index of tng crv model: NrlMdlTrpl[3]
        varTmpX = int(np.around(NrlMdlTrpl[0], 0))
        varTmpY = int(np.around(NrlMdlTrpl[1], 0))

        # Create pRF model (2D):
        aryGauss = funcGauss2D(varPixX,
                               varPixY,
                               varTmpX,
                               varTmpY,
                               NrlMdlTrpl[2])

        # Multiply pixel-wise box car model with Gaussian pRF models:
        aryNrlTcTmp = np.multiply(aryBoxCar, aryGauss[:, :, None, None])

        # Calculate sum across x- and y-dimensions - the 'area under the
        # Gaussian surface'. This is essentially an unscaled version of the
        # neural time course model (i.e. not yet scaled for the size of
        # the pRF).
        aryNrlTcTmp = np.sum(aryNrlTcTmp, axis=(0, 1))

        # Normalise the nrl time course model to the size of the pRF. This
        # gives us the ratio of 'activation' of the pRF at each time point,
        # or, in other words, the neural time course model.
        aryNrlTcTmp = np.divide(aryNrlTcTmp,
                                np.sum(aryGauss, axis=(0, 1)))

        # Put model time courses into the function's output array:
        aryOut[idx, :, :] = aryNrlTcTmp

        # Status indicator (only used in the first of the parallel
        # processes):
        if idxPrc == 1:
            # Increment status indicator counter:
            varCntSts02 = varCntSts02 + 1

    # Output list:
    lstOut = [idxPrc,
              aryOut,
              ]

    queOut.put(lstOut)
841778
def process_mav(self, mlog, flightmode_selections):
        ''''''
        self.vars = {}
        idx = 0
        all_false = True
        for s in flightmode_selections:
            if s:
                all_false = False

        # pre-calc right/left axes
        self.num_fields = len(self.fields)
        for i in range(0, self.num_fields):
            f = self.fields[i]
            if f.endswith(":2"):
                self.axes[i] = 2
                f = f[:-2]
            if f.endswith(":1"):
                self.first_only[i] = True
                f = f[:-2]
            self.fields[i] = f

        # see which fields are simple
        self.simple_field = []
        for i in range(0, self.num_fields):
            f = self.fields[i]
            m = re.match('^([A-Z][A-Z0-9_]*)[.]([A-Za-z_][A-Za-z0-9_]*)$', f)
            if m is None:
                self.simple_field.append(None)
            else:
                self.simple_field.append((m.group(1),m.group(2)))

        if len(self.flightmode_list) > 0:
            # prime the timestamp conversion
            self.timestamp_to_days(self.flightmode_list[0][1])

        while True:
            msg = mlog.recv_match(type=self.msg_types)
            if msg is None:
                break
            if msg.get_type() not in self.msg_types:
                continue
            if self.condition:
                if not mavutil.evaluate_condition(self.condition, mlog.messages):
                    continue
            tdays = self.timestamp_to_days(msg._timestamp)

            if all_false or len(flightmode_selections) == 0:
                self.add_data(tdays, msg, mlog.messages)
            else:
                if idx < len(self.flightmode_list) and msg._timestamp >= self.flightmode_list[idx][2]:
                    idx += 1
                elif (idx < len(flightmode_selections) and flightmode_selections[idx]):
                    self.add_data(tdays, msg, mlog.messages)
1363014
def deltas(last, now):
    """"""
    return {
        xy: RouterDiagnostics(*((n - l) & 0xFFFFFFFF
                                for l, n in zip(last[xy], now[xy])))
        for xy in last
    }
301871
private void estimateJitter(RtpPacket newPacket) {
        long arrival = rtpClock.getLocalRtpTime();
        long newPacketTimestamp = newPacket.getTimestamp();
        long transit = arrival - newPacketTimestamp;
        long d = transit - currentTransit;
        if (d < 0) {
            d = -d;
        }
        // logger.info(String.format("recalculating jitter: arrival=%d, newPacketTimestamp=%d, transit=%d, transit delta=%d",
        // arrival, newPacketTimestamp, transit, d ));
        currentTransit = transit;
        currentJitter += d - ((currentJitter + 8) >> 4);
    }
127092
def k_modes(X, n_clusters, max_iter, dissim, init, n_init, verbose, random_state, n_jobs):
    """"""
    random_state = check_random_state(random_state)
    if sparse.issparse(X):
        raise TypeError("k-modes does not support sparse data.")

    X = check_array(X, dtype=None)

    # Convert the categorical values in X to integers for speed.
    # Based on the unique values in X, we can make a mapping to achieve this.
    X, enc_map = encode_features(X)

    n_points, n_attrs = X.shape
    assert n_clusters <= n_points, "Cannot have more clusters ({}) " \
                                   "than data points ({}).".format(n_clusters, n_points)

    # Are there more n_clusters than unique rows? Then set the unique
    # rows as initial values and skip iteration.
    unique = get_unique_rows(X)
    n_unique = unique.shape[0]
    if n_unique <= n_clusters:
        max_iter = 0
        n_init = 1
        n_clusters = n_unique
        init = unique

    results = []
    seeds = random_state.randint(np.iinfo(np.int32).max, size=n_init)
    if n_jobs == 1:
        for init_no in range(n_init):
            results.append(k_modes_single(X, n_clusters, n_points, n_attrs, max_iter,
                                          dissim, init, init_no, verbose, seeds[init_no]))
    else:
        results = Parallel(n_jobs=n_jobs, verbose=0)(
            delayed(k_modes_single)(X, n_clusters, n_points, n_attrs, max_iter,
                                    dissim, init, init_no, verbose, seed)
            for init_no, seed in enumerate(seeds))
    all_centroids, all_labels, all_costs, all_n_iters = zip(*results)

    best = np.argmin(all_costs)
    if n_init > 1 and verbose:
        print("Best run was number {}".format(best + 1))

    return all_centroids[best], enc_map, all_labels[best], \
        all_costs[best], all_n_iters[best]
1308399
def doublewell_eigs(n_grid, lag_time=1):
    
    return _brownian_eigs(n_grid, lag_time, DOUBLEWELL_GRAD_POTENTIAL,
                          -np.pi, np.pi, reflect_bc=True)