Attach timeseries to indicators, for performing Bayesian inference.
338825
def get_timeseries_values_for_indicators(
        self, resolution: str = "month", months: Iterable[int] = range(6, 9)
    ):
        
        if resolution == "month":
            funcs = [
                partial(get_indicator_value, month=month) for month in months
            ]
        else:
            raise NotImplementedError(
                "Currently, only the 'month' resolution is supported."
            )

        for n in self.nodes(data=True):
            for indicator in n[1]["indicators"].values():
                indicator.timeseries = [
                    func(indicator, year="2017")[0] for func in funcs
                ]
                if len(set(indicator.timeseries)) == 1:
                    indicator.timeseries = None
1306392
def make_ica_funs(observed_dimension, latent_dimension):
    """"""

    def sample(weights, n_samples, noise_std, rs):
        latents = rs.randn(latent_dimension, n_samples)
        latents = np.array(sorted(latents.T, key=lambda a_entry: a_entry[0])).T
        noise = rs.randn(n_samples, observed_dimension) * noise_std
        observed = predict(weights, latents) + noise
        return latents, observed

    def predict(weights, latents):
        return np.dot(weights, latents).T

    def logprob(weights, latents, noise_std, observed):
        preds = predict(weights, latents)
        log_lik = np.sum(t.logpdf(preds, 2.4, observed, noise_std))
        return log_lik

    num_weights = observed_dimension * latent_dimension

    def unpack_weights(weights):
        return np.reshape(weights, (observed_dimension, latent_dimension))

    return num_weights, sample, logprob, unpack_weights
917828
def inference_minibatch(self, kern, X, Z, likelihood, Y, Zp):
        

        num_data, output_dim = Y.shape

        #see whether we've got a different noise variance for each datum
        beta = 1./np.fmax(likelihood.variance, 1e-6)
        het_noise = beta.size > 1
        # VVT_factor is a matrix such that tdot(VVT_factor) = VVT...this is for efficiency!
        #self.YYTfactor = beta*self.get_YYTfactor(Y)
        if self.Y_speedup and not het_noise:
            YYT_factor = self.get_YYTfactor(Y)
        else:
            YYT_factor = Y

        n_start = self.batch_pos
        batchsize = num_data if self.batchsize is None else self.batchsize
        n_end = min(batchsize+n_start, num_data)
        if n_end == num_data:
            isEnd = True
            self.batch_pos = 0
        else:
            isEnd = False
            self.batch_pos = n_end

        if batchsize == num_data:
            Y_slice = YYT_factor
            X_slice = X
        else:
            Y_slice = YYT_factor[n_start:n_end]
            X_slice = X[n_start:n_end]

        psi0 = kern._Kdiag(X_slice) #Kffdiag
        psi1 = kern.K(X_slice, Z) #Kfu
        betapsi1 = np.einsum('n,nm->nm', beta, psi1)

        X_slice = X_slice.values
        Z = Z.values

        Zp = Zp.gamma.values
        indX = np.int_(X_slice[:, -1])
        indZ = np.int_(Z[:, -1]) - Zp.shape[0]

        betaY = beta*Y_slice

        #======================================================================
        # Load Intermediate Results
        #======================================================================

        dL_dpsi2R = self.midRes['dL_dpsi2R']
        v = self.midRes['v']

        #======================================================================
        # Compute dL_dpsi
        #======================================================================

        dL_dpsi0 = -.5*output_dim*(beta * Zp[indX, :]) #XxQ #TODO: Check this gradient

        dL_dpsi1 = np.dot(betaY, v.T)
        dL_dEZp  = psi1*dL_dpsi1
        dL_dpsi1 = Zp[np.ix_(indX, indZ)]*dL_dpsi1
        dL_dgamma = np.zeros(Zp.shape)
        for d in np.unique(indX):
            indd = indX == d
            betapsi1d = betapsi1[indd, :]
            psi1d = psi1[indd, :]
            Zpd = Zp[d, :]
            Zp2 = Zpd[:, None]*Zpd[None, :] - np.diag(np.power(Zpd, 2)) + np.diag(Zpd)
            dL_dpsi1[indd, :] += np.dot(betapsi1d, Zp2[np.ix_(indZ, indZ)] * dL_dpsi2R)*2.

            dL_EZp2 = dL_dpsi2R * (np.dot(psi1d.T, psi1d) * beta)*2.  # Zpd*Kufd*Kfud*beta
            #Gradient of Likelihood wrt gamma is calculated here
            EZ = Zp[d, indZ]
            for q in range(Zp.shape[1]):
                EZt = EZ.copy()
                indq = indZ == q
                EZt[indq] = .5
                dL_dgamma[d, q] = np.sum(dL_dEZp[np.ix_(indd, indq)]) + np.sum(dL_EZp2[:, indq]*EZt[:, None]) -\
                    .5*beta*(np.sum(psi0[indd, q]))

        #======================================================================
        # Compute dL_dthetaL
        #======================================================================
        if isEnd:
            dL_dthetaL = self.midRes['dL_dthetaL']
        else:
            dL_dthetaL = 0.

        grad_dict = {'dL_dKdiag': dL_dpsi0,
                     'dL_dKnm': dL_dpsi1,
                     'dL_dthetaL': dL_dthetaL,
                     'dL_dgamma': dL_dgamma}

        return isEnd, (n_start, n_end), grad_dict
605274
def _construct_timeseries(self, timeseries, constraints={}):
        
        self.response_from(timeseries, constraints)
        if self.response == None:
            return None
        return {'data':self.response['data'],
                'period':self.response['period'],
                'start time':datetime.datetime.fromtimestamp(self.response['start_time']),
                'end time':datetime.datetime.fromtimestamp(self.response['end_time'])}
86810
def generate_values(self, *args, **kwargs):
        
        assert 'ref value' in self.kwargs

        # 1. Generate $\mu$
        start_date = self.times[0].to_pydatetime()
        end_date = self.times[-1].to_pydatetime()
        ref_date = self.ref_date

        mu = self.generate_mu(end_date, ref_date, start_date)

        # 3. Generate $\sigma$
        ## Prepare array with growth values $\sigma$
        if self.sample_mean_value:
            sigma = np.zeros((len(self.times), self.size))
        else:

            if self.kwargs['type'] == 'interp':

                def get_date(record):
                    return datetime.datetime.strptime(record[0], "%Y-%m-%d")

                ref_value_ = sorted(json.loads(self.kwargs['ref value'].strip()).items(), key=get_date)
                intial_value = ref_value_[0][1]
            else:
                intial_value = self.kwargs['ref value']

            variability_ = intial_value * self.kwargs['initial_value_proportional_variation']
            logger.debug(f'sampling random distribution with parameters -{variability_}, 0, {variability_}')
            sigma = np.random.triangular(-1 * variability_, 0, variability_,
                                         (len(self.times), self.size))
        ## 4. Prepare growth array for $\alpha_{sigma}$
        alpha_sigma = growth_coefficients(start_date,
                                          end_date,
                                          ref_date,
                                          self.kwargs['ef_growth_factor'], 1)

        ### 5. Prepare DataFrame
        iterables = [self.times, range(self.size)]
        index_names = ['time', 'samples']
        _multi_index = pd.MultiIndex.from_product(iterables, names=index_names)

        df = pd.DataFrame(index=_multi_index, dtype=float)

        from dateutil import relativedelta
        r = relativedelta.relativedelta(end_date, start_date)
        months = r.years * 12 + r.months + 1
        name = kwargs['name']
        ## Apply growth to $\sigma$ and add $\sigma$ to $\mu$
        df[name] = ((sigma * alpha_sigma) + mu.reshape(months, 1)).ravel()

        ## test if df has sub-zero values
        df_sigma__dropna = df[name].where(df[name] < 0).dropna()
        if not df_sigma__dropna.empty:
            logger.warning(f"Negative values for parameter {name} from {df_sigma__dropna.index[0][0]}")

        return df[name]
441877
def util_granulate_time_series(time_series, scale):
    
    n = len(time_series)
    b = int(np.fix(n / scale))
    temp = np.reshape(time_series[0:b*scale], (b, scale))
    cts = np.mean(temp, axis = 1)
    return cts
63479
def event_estimation(self, event, logliks, logsumexps, renamed=''):
        
        plotter = _ModelLoglikPlotter()
        plotter.plot(event, logliks, logsumexps, self.modality_to_color,
                     renamed=renamed)
        return plotter
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
867655
private List<SequenceState<State, Observation, Path>> computeViterbiSequence(
            List<TimeStep<State, Observation, Path>> timeSteps, int originalGpxEntriesCount,
            QueryGraph queryGraph) {
        final HmmProbabilities probabilities
                = new HmmProbabilities(measurementErrorSigma, transitionProbabilityBeta);
        final ViterbiAlgorithm<State, Observation, Path> viterbi = new ViterbiAlgorithm<>();

        logger.debug("\n=============== Paths ===============");
        int timeStepCounter = 0;
        TimeStep<State, Observation, Path> prevTimeStep = null;
        int i = 1;
        for (TimeStep<State, Observation, Path> timeStep : timeSteps) {
            logger.debug("\nPaths to time step {}", i++);
            computeEmissionProbabilities(timeStep, probabilities);

            if (prevTimeStep == null) {
                viterbi.startWithInitialObservation(timeStep.observation, timeStep.candidates,
                        timeStep.emissionLogProbabilities);
            } else {
                computeTransitionProbabilities(prevTimeStep, timeStep, probabilities, queryGraph);
                viterbi.nextStep(timeStep.observation, timeStep.candidates,
                        timeStep.emissionLogProbabilities, timeStep.transitionLogProbabilities,
                        timeStep.roadPaths);
            }
            if (viterbi.isBroken()) {
                String likelyReasonStr = "";
                if (prevTimeStep != null) {
                    Observation prevGPXE = prevTimeStep.observation;
                    Observation gpxe = timeStep.observation;
                    double dist = distanceCalc.calcDist(prevGPXE.getPoint().lat, prevGPXE.getPoint().lon,
                            gpxe.getPoint().lat, gpxe.getPoint().lon);
                    if (dist > 2000) {
                        likelyReasonStr = "Too long distance to previous measurement? "
                                + Math.round(dist) + "m, ";
                    }
                }

                throw new IllegalArgumentException("Sequence is broken for submitted track at time step "
                        + timeStepCounter + " (" + originalGpxEntriesCount + " points). "
                        + likelyReasonStr + "observation:" + timeStep.observation + ", "
                        + timeStep.candidates.size() + " candidates: "
                        + getSnappedCandidates(timeStep.candidates)
                        + ". If a match is expected consider increasing max_visited_nodes.");
            }

            timeStepCounter++;
            prevTimeStep = timeStep;
        }

        return viterbi.computeMostLikelySequence();
    }
354029
def _plot_posterior_op(
    values,
    var_name,
    selection,
    ax,
    bw,
    linewidth,
    bins,
    kind,
    point_estimate,
    round_to,
    credible_interval,
    ref_val,
    rope,
    ax_labelsize,
    xt_labelsize,
    **kwargs
):  # noqa: D202
    """"""

    def format_as_percent(x, round_to=0):
        return "{0:.{1:d}f}%".format(100 * x, round_to)

    def display_ref_val():
        if ref_val is None:
            return
        elif isinstance(ref_val, dict):
            val = None
            for sel in ref_val.get(var_name, []):
                if all(
                    k in selection and selection[k] == v for k, v in sel.items() if k != "ref_val"
                ):
                    val = sel["ref_val"]
                    break
            if val is None:
                return
        elif isinstance(ref_val, Number):
            val = ref_val
        else:
            raise ValueError(
                "Argument `ref_val` must be None, a constant, or a "
                'dictionary like {"var_name": {"ref_val": (lo, hi)}}'
            )
        less_than_ref_probability = (values < val).mean()
        greater_than_ref_probability = (values >= val).mean()
        ref_in_posterior = "{} <{:g}< {}".format(
            format_as_percent(less_than_ref_probability, 1),
            val,
            format_as_percent(greater_than_ref_probability, 1),
        )
        ax.axvline(val, ymin=0.05, ymax=0.75, color="C1", lw=linewidth, alpha=0.65)
        ax.text(
            values.mean(),
            plot_height * 0.6,
            ref_in_posterior,
            size=ax_labelsize,
            color="C1",
            weight="semibold",
            horizontalalignment="center",
        )

    def display_rope():
        if rope is None:
            return
        elif isinstance(rope, dict):
            vals = None
            for sel in rope.get(var_name, []):
                if all(k in selection and selection[k] == v for k, v in sel.items() if k != "rope"):
                    vals = sel["rope"]
                    break
            if vals is None:
                return
        elif len(rope) == 2:
            vals = rope
        else:
            raise ValueError(
                "Argument `rope` must be None, a dictionary like"
                '{"var_name": {"rope": (lo, hi)}}, or an'
                "iterable of length 2"
            )

        ax.plot(
            vals,
            (plot_height * 0.02, plot_height * 0.02),
            lw=linewidth * 5,
            color="C2",
            solid_capstyle="round",
            zorder=0,
            alpha=0.7,
        )
        text_props = {"size": ax_labelsize, "horizontalalignment": "center", "color": "C2"}
        ax.text(vals[0], plot_height * 0.2, vals[0], weight="semibold", **text_props)
        ax.text(vals[1], plot_height * 0.2, vals[1], weight="semibold", **text_props)

    def display_point_estimate():
        if not point_estimate:
            return
        if point_estimate not in ("mode", "mean", "median"):
            raise ValueError("Point Estimate should be in ('mode','mean','median')")
        if point_estimate == "mean":
            point_value = values.mean()
        elif point_estimate == "mode":
            if isinstance(values[0], float):
                density, lower, upper = _fast_kde(values, bw=bw)
                x = np.linspace(lower, upper, len(density))
                point_value = x[np.argmax(density)]
            else:
                point_value = mode(values.round(round_to))[0][0]
        elif point_estimate == "median":
            point_value = np.median(values)
        point_text = "{}={:.{}f}".format(point_estimate, point_value, round_to)

        ax.text(
            point_value,
            plot_height * 0.8,
            point_text,
            size=ax_labelsize,
            horizontalalignment="center",
        )

    def display_hpd():
        hpd_intervals = hpd(values, credible_interval=credible_interval)
        ax.plot(
            hpd_intervals,
            (plot_height * 0.02, plot_height * 0.02),
            lw=linewidth * 2,
            color="k",
            solid_capstyle="round",
        )
        ax.text(
            hpd_intervals[0],
            plot_height * 0.07,
            hpd_intervals[0].round(round_to),
            size=ax_labelsize,
            horizontalalignment="center",
        )
        ax.text(
            hpd_intervals[1],
            plot_height * 0.07,
            hpd_intervals[1].round(round_to),
            size=ax_labelsize,
            horizontalalignment="center",
        )
        ax.text(
            (hpd_intervals[0] + hpd_intervals[1]) / 2,
            plot_height * 0.3,
            format_as_percent(credible_interval) + " HPD",
            size=ax_labelsize,
            horizontalalignment="center",
        )

    def format_axes():
        ax.yaxis.set_ticks([])
        ax.spines["top"].set_visible(False)
        ax.spines["right"].set_visible(False)
        ax.spines["left"].set_visible(False)
        ax.spines["bottom"].set_visible(True)
        ax.xaxis.set_ticks_position("bottom")
        ax.tick_params(
            axis="x", direction="out", width=1, length=3, color="0.5", labelsize=xt_labelsize
        )
        ax.spines["bottom"].set_color("0.5")

    if kind == "kde" and values.dtype.kind == "f":
        plot_kde(
            values,
            bw=bw,
            fill_kwargs={"alpha": kwargs.pop("fill_alpha", 0)},
            plot_kwargs={"linewidth": linewidth},
            ax=ax,
            rug=False,
        )
    else:
        if bins is None:
            if values.dtype.kind == "i":
                xmin = values.min()
                xmax = values.max()
                bins = range(xmin, xmax + 2)
                ax.set_xlim(xmin - 0.5, xmax + 0.5)
            else:
                bins = "auto"
        kwargs.setdefault("align", "left")
        kwargs.setdefault("color", "C0")
        ax.hist(values, bins=bins, alpha=0.35, **kwargs)

    plot_height = ax.get_ylim()[1]

    format_axes()
    display_hpd()
    display_point_estimate()
    display_ref_val()
    display_rope()
746226
def stick_bgplvm(model=None, optimize=True, verbose=True, plot=True):
    """"""
    from GPy.models import BayesianGPLVM
    from matplotlib import pyplot as plt
    import numpy as np
    import GPy
    import pods

    data = pods.datasets.osu_run1()
    Q = 6
    kernel = GPy.kern.RBF(Q, lengthscale=np.repeat(.5, Q), ARD=True)
    m = BayesianGPLVM(data['Y'], Q, init="PCA", num_inducing=20, kernel=kernel)

    m.data = data
    m.likelihood.variance = 0.001

    # optimize
    try:
        if optimize: m.optimize('bfgs', messages=verbose, max_iters=5e3, bfgs_factor=10)
    except KeyboardInterrupt:
        print("Keyboard interrupt, continuing to plot and return")

    if plot:
        fig, (latent_axes, sense_axes) = plt.subplots(1, 2)
        plt.sca(latent_axes)
        m.plot_latent(ax=latent_axes)
        y = m.Y[:1, :].copy()
        data_show = GPy.plotting.matplot_dep.visualize.stick_show(y, connect=data['connect'])
        dim_select = GPy.plotting.matplot_dep.visualize.lvm_dimselect(m.X.mean[:1, :].copy(), m, data_show, latent_axes=latent_axes, sense_axes=sense_axes)
        fig.canvas.draw()
        # Canvas.show doesn't work on OSX.
        #fig.canvas.show()
        input('Press enter to finish')

    return m
227310
def estimate_map(interface, state, label, inp):
    """"""
    out = interface.output(0)
    centers = {}
    for row in inp:
        row = row.strip().split(state["delimiter"])
        if len(row) > 1:
            x = [(0 if row[i] in state["missing_vals"] else float(row[i])) for i in state["X_indices"]]
            cluster = min((state['dist'](c, x), i) for i, c in state['centers'])[1]
            vertex = state['create'](x, 1.0)
            centers[cluster] = vertex if cluster not in centers else state["update"](centers[cluster], vertex)

    for cluster, values in centers.iteritems():
        out.add(cluster, values)
202025
def _update_inernal(self, varp_list):
        """"""
        # The probability for the binary variable for the same latent dimension of any of the models is on.
        if self.group_spike:
            self._b_prob_all = 1.-param_to_array(varp_list[0].gamma_group)
            [np.multiply(self._b_prob_all, 1.-vp.gamma_group, self._b_prob_all) for vp in varp_list[1:]]
        else:
            self._b_prob_all = 1.-param_to_array(varp_list[0].binary_prob)
            [np.multiply(self._b_prob_all, 1.-vp.binary_prob, self._b_prob_all) for vp in varp_list[1:]]
944469
def cumulative_density_at_times(self, times, label=None):
        
        label = coalesce(label, self._label)
        return pd.Series(1 - self.predict(times), index=_to_array(times), name=label)
579729
def _check_time_series(self, other):
        
        if not isinstance(other, TimeSeries):
            msg = "unsupported operand types(s) for +: %s and %s" % \
                  (type(self), type(other))
            raise TypeError(msg)
656272
def compute(cls, observation, prediction):
        """"""
        assert isinstance(observation, dict)
        try:
            p_value = prediction['mean']  # Use the prediction's mean.
        except (TypeError, KeyError, IndexError):  # If there isn't one...
            try:
                p_value = prediction['value']  # Use the prediction's value.
            except (TypeError, IndexError):  # If there isn't one...
                p_value = prediction  # Use the prediction (assume numeric).
        o_mean = observation['mean']
        o_std = observation['std']
        value = (p_value - o_mean)/o_std
        value = utils.assert_dimensionless(value)
        if np.isnan(value):
            score = InsufficientDataScore('One of the input values was NaN')
        else:
            score = ZScore(value)
        return score
1069195
def inference(self, observed_arr):
        '''
        
        '''
        _ = self.__lstm_model.inference(observed_arr)
        return self.__lstm_model.get_feature_points()