Apply correlation to randomly sampled residuals.

        See Parent function
8049
def apply_correlation(self, sites, imt, residuals, stddev_intra):
        
        # stddev_intra is repeated if it is only 1 value for all the residuals
        if stddev_intra.shape[0] == 1:
            stddev_intra = numpy.matlib.repmat(
                stddev_intra, len(sites.complete), 1)
        # Reshape 'stddev_intra' if needed
        stddev_intra = stddev_intra.squeeze()
        if not stddev_intra.shape:
            stddev_intra = stddev_intra[None]

        if self.uncertainty_multiplier == 0:   # No uncertainty

            # residuals were sampled from a normal distribution with
            # stddev_intra standard deviation. 'residuals_norm' are residuals
            # normalized, sampled from a standard normal distribution.
            # For this, every row of 'residuals' (every site) is divided by its
            # corresponding standard deviation element.
            residuals_norm = residuals / stddev_intra[sites.sids, None]

            # Lower diagonal of the Cholesky decomposition from/to cache
            try:
                cormaLow = self.cache[imt]
            except KeyError:
                # Note that instead of computing the whole correlation matrix
                # corresponding to sites.complete, here we compute only the
                # correlation matrix corresponding to sites.
                cormaLow = numpy.linalg.cholesky(
                       numpy.diag(stddev_intra[sites.sids]) *
                       self._get_correlation_matrix(sites, imt) *
                       numpy.diag(stddev_intra[sites.sids]))
                self.cache[imt] = cormaLow

            # Apply correlation
            return numpy.dot(cormaLow, residuals_norm)

        else:   # Variability (uncertainty) is included
            nsim = len(residuals[1])
            nsites = len(residuals)

            # Re-sample all the residuals
            residuals_correlated = residuals * 0
            for isim in range(0, nsim):
                corma = self._get_correlation_matrix(sites, imt)
                cov = (numpy.diag(stddev_intra[sites.sids]) * corma *
                       numpy.diag(stddev_intra[sites.sids]))
                residuals_correlated[0:, isim] = (
                    numpy.random.multivariate_normal(
                        numpy.zeros(nsites), cov, 1))

            return residuals_correlated
327156
def _refit(self):
        
        # temporary variables for manipulation
        mu_diff = np.empty_like(self.f.components[0].mu)
        sigma   = np.empty_like(self.f.components[0].sigma)
        mean    = np.empty_like(mu_diff)
        cov     = np.empty_like(sigma)

        for j, c in enumerate(self.g.components):
            # stop if inv_map is empty for j-th comp.
            if not self.inv_map[j]:
                self.g.weights[j] = 0.
                continue

            # (re-)initialize new mean/cov to zero
            mean[:] = 0.0
            cov[:] = 0.0

            # compute total weight and mean
            self.g.weights[j] = self.f.weights[self.inv_map[j]].sum()
            for i in self.inv_map[j]:
                mean += self.f.weights[i] * self.f.components[i].mu

            # rescale by total weight
            mean /= self.g.weights[j]

            # update covariance
            for i in self.inv_map[j]:
                # mu_diff = mu'_j - mu_i
                mu_diff[:] = mean
                mu_diff -= self.f.components[i].mu

                # sigma = (mu'_j - mu_i) (mu'_j - mu_i)^T
                sigma[:] = np.outer(mu_diff, mu_diff)

                # sigma += sigma_i
                sigma += self.f.components[i].sigma

                # multiply with alpha_i
                sigma *= self.f.weights[i]

                # sigma_j += alpha_i * (sigma_i + (mu'_j - mu_i) (mu'_j - mu_i)^T
                cov += sigma

            # 1 / beta_j
            cov /= self.g.weights[j]

            # update the Mixture
            c.update(mean, cov)
630441
def generate_covariance(ts, sigma, tau):

    
    """

    ndim = ts.shape[0]

    tis = ts[:, np.newaxis]
    tjs = ts[np.newaxis, :]

    return sigma*sigma*np.exp(-np.square(tis-tjs)/(2.0*tau*tau))
1038395
def _sample_row(self):
        
        unis = np.random.uniform(0, 1, self.n_var)
        # randomly select a node to start with
        first_ind = np.random.randint(0, self.n_var)
        adj = self.trees[0].get_adjacent_matrix()
        visited = []
        explore = [first_ind]

        sampled = np.zeros(self.n_var)
        itr = 0
        while explore:
            current = explore.pop(0)
            neighbors = np.where(adj[current, :] == 1)[0].tolist()
            if itr == 0:
                new_x = self.ppfs[current](unis[current])

            else:
                for i in range(itr - 1, -1, -1):
                    current_ind = -1

                    if i >= self.truncated:
                        continue

                    current_tree = self.trees[i].edges
                    # get index of edge to retrieve
                    for edge in current_tree:
                        if i == 0:
                            if (edge.L == current and edge.R == visited[0]) or\
                               (edge.R == current and edge.L == visited[0]):
                                current_ind = edge.index
                                break
                        else:
                            if edge.L == current or edge.R == current:
                                condition = set(edge.D)
                                condition.add(edge.L)
                                condition.add(edge.R)

                                visit_set = set(visited)
                                visit_set.add(current)

                                if condition.issubset(visit_set):
                                    current_ind = edge.index
                                break

                    if current_ind != -1:
                        # the node is not indepedent contional on visited node
                        copula_type = current_tree[current_ind].name
                        copula = Bivariate(CopulaTypes(copula_type))
                        copula.theta = current_tree[current_ind].theta
                        derivative = copula.partial_derivative_scalar

                        if i == itr - 1:
                            tmp = optimize.fminbound(
                                derivative, EPSILON, 1.0,
                                args=(unis[visited[0]], unis[current])
                            )
                        else:
                            tmp = optimize.fminbound(
                                derivative, EPSILON, 1.0,
                                args=(unis[visited[0]], tmp)
                            )

                        tmp = min(max(tmp, EPSILON), 0.99)

                new_x = self.ppfs[current](tmp)

            sampled[current] = new_x

            for s in neighbors:
                if s not in visited:
                    explore.insert(0, s)

            itr += 1
            visited.insert(0, current)

        return sampled
914991
def resample(self, destination=None, **kwargs):
        """"""
        return self._generate_scene_func(self._scenes, 'resample', True, destination=destination, **kwargs)
413876
def plot_series_residuals_internal(self, varied_data, varied_idx, **kwargs):
        
        nf = len(self.f_cb(*self.pre_process(
            self.internal_xout[0], self.internal_params_out[0])))
        xerr = np.empty((self.internal_xout.shape[0], nf))
        for idx, (res, params) in enumerate(zip(self.internal_xout, self.internal_params_out)):
            xerr[idx, :] = self.f_cb(res, params)
        return self.plot_series(xerr, varied_data, varied_idx, **kwargs)
327681
def keep_negative_resample(X, y, model_generator, method_name, num_fcounts=11):
    
    return __run_measure(measures.keep_resample, X, y, model_generator, method_name, -1, num_fcounts, __mean_pred)
1292126
def par_compute_residuals(i):
    

    # Compute the residuals in parallel, need to check if the residuals
    # depend on alpha
    global mp_ry0
    global mp_ry1
    global mp_sy0
    global mp_sy1
    global mp_nrmAx
    global mp_nrmBy
    global mp_nrmu
    mp_ry0[i] = np.sum((mp_DXnr[i] - mp_Y0[i])**2)
    mp_ry1[i] = mp_alpha**2*np.sum((mp_Xnr[mp_grp[i]:mp_grp[i+1]]-
                                    mp_Y1[mp_grp[i]:mp_grp[i+1]])**2)
    mp_sy0[i] = np.sum((mp_Y0old[i] - mp_Y0[i])**2)
    mp_sy1[i] = mp_alpha**2*np.sum((mp_Y1old[mp_grp[i]:mp_grp[i+1]]-
                                    mp_Y1[mp_grp[i]:mp_grp[i+1]])**2)
    mp_nrmAx[i] = np.sum(mp_DXnr[i]**2) + mp_alpha**2 * np.sum(
        mp_Xnr[mp_grp[i]:mp_grp[i+1]]**2)
    mp_nrmBy[i] = np.sum(mp_Y0[i]**2) + mp_alpha**2 * np.sum(
        mp_Y1[mp_grp[i]:mp_grp[i+1]]**2)
    mp_nrmu[i] = np.sum(mp_U0[i]**2) + np.sum(mp_U1[mp_grp[i]:mp_grp[i+1]]**2)
565720
def resampling_scheme(func):
    """"""

    @functools.wraps(func)
    def modif_func(W, M=None):
        M = W.shape[0] if M is None else M
        return func(W, M)

    rs_funcs[func.__name__] = modif_func
    modif_func.__doc__ = rs_doc % func.__name__.capitalize()
    return modif_func
993443
def Residual(*layers, **kwargs):
  """"""
  shortcut = kwargs.get('shortcut', Identity())  # pylint: disable=no-value-for-parameter
  if len(layers) > 1:
    return Serial(
        Branch(),  # pylint: disable=no-value-for-parameter
        Parallel(Serial(*layers), shortcut),
        SumBranches()  # pylint: disable=no-value-for-parameter
    )
  elif len(layers) == 1:
    return Serial(
        Branch(),  # pylint: disable=no-value-for-parameter
        Parallel(layers[0], shortcut),
        SumBranches()  # pylint: disable=no-value-for-parameter
    )
  else:
    raise ValueError('Empty residual combinator.')
1167176
def cv_residuals(self, cv=True):
        """"""
        vals = self.cv_values(cv)
        return (self.y - vals) / self.dy
200404
private static double[] calcLinReg(double[] xData, double[] yData) {
        double sumX = 0;
        double sumY = 0;
        double prod = 0;
        int nstat = xData.length;
        double[] regCoef = new double[3]; //(intercept, gradient, r2)

        double meanYValue = mean(yData);
        double meanXValue = mean(xData);

        //calculating regression coefficients
        for (int i = 0; i < nstat; i++) {
            sumX += (xData[i] - meanXValue) * (xData[i] - meanXValue);
            sumY += (yData[i] - meanYValue) * (yData[i] - meanYValue);
            prod += (xData[i] - meanXValue) * (yData[i] - meanYValue);
        }
        if (sumX > 0 && sumY > 0) {
            regCoef[1] = prod / sumX;  //gradient
            regCoef[0] = meanYValue - regCoef[1] * meanXValue; //intercept
            regCoef[2] = Math.pow((prod / Math.sqrt(sumX * sumY)), 2); //r2
        }
        return regCoef;
    }
197318
def resample(self, data, mask=None, omega=None):
        
        if not isinstance(data, list):
            assert isinstance(data, tuple) and len(data) == 2, \
                "datas must be an (x,y) tuple or a list of such tuples"
            data = [data]

        if mask is None:
            mask = [np.ones(y.shape, dtype=bool) for x, y in data]

        # Resample auxiliary variables if they are not given
        if omega is None:
            omega = self._resample_auxiliary_variables(data)

        # Make copies of parameters (for sample collection in calling methods)
        self.A = self.A.copy()
        self.b = self.b.copy()

        D = self.D_in
        for n in range(self.D_out):
            # Resample C_{n,:} given z, omega[:,n], and kappa[:,n]
            prior_Sigma = np.zeros((D + 1, D + 1))
            prior_Sigma[:D, :D] = self.sigmasq_A[n]
            prior_Sigma[D, D] = self.sigmasq_b[n]
            prior_J = np.linalg.inv(prior_Sigma)

            prior_h = prior_J.dot(np.concatenate((self.mu_A[n], [self.mu_b[n]])))

            lkhd_h = np.zeros(D + 1)
            lkhd_J = np.zeros((D + 1, D + 1))

            for d, m, o in zip(data, mask, omega):
                if isinstance(d, tuple):
                    x, y = d
                else:
                    x, y = d[:, :D], d[:, D:]
                augx = np.hstack((x, np.ones((x.shape[0], 1))))
                J = o * m
                h = self.kappa_func(y) * m

                lkhd_J += (augx * J[:, n][:, None]).T.dot(augx)
                lkhd_h += h[:, n].T.dot(augx)

            post_h = prior_h + lkhd_h
            post_J = prior_J + lkhd_J

            joint_sample = sample_gaussian(J=post_J, h=post_h)
            self.A[n, :] = joint_sample[:D]
            self.b[n] = joint_sample[D]
910304
def compute_residuals(self):
        

        # The full primary residual is straightforward to compute from
        # the primary residuals for the baseline algorithm and for the
        # additional variables
        r0 = self.rsdl_r(self.AXnr, self.Y)
        r1 = self.AX1nr - self.Y1 - self.S
        r = np.sqrt(np.sum(r0**2) + np.sum(r1**2))

        # The full dual residual is more complicated to compute than the
        # full primary residual
        ATU = self.swapaxes(self.U) + sl.irfftn(
            np.conj(self.Zf) * sl.rfftn(self.U1, self.cri.Nv, self.cri.axisN),
            self.cri.Nv, self.cri.axisN)
        s = self.rho * np.linalg.norm(ATU)

        # The normalisation factor for the full primal residual is also not
        # straightforward
        nAX = np.sqrt(np.linalg.norm(self.AXnr)**2 +
                      np.linalg.norm(self.AX1nr)**2)
        nY = np.sqrt(np.linalg.norm(self.Y)**2 +
                     np.linalg.norm(self.Y1)**2)
        rn = max(nAX, nY, np.linalg.norm(self.S))

        # The normalisation factor for the full dual residual is
        # straightforward to compute
        sn = self.rho * np.sqrt(np.linalg.norm(self.U)**2 +
                                np.linalg.norm(self.U1)**2)

        # Final residual values and stopping tolerances depend on
        # whether standard or normalised residuals are specified via the
        # options object
        if self.opt['AutoRho', 'StdResiduals']:
            epri = np.sqrt(self.Nc)*self.opt['AbsStopTol'] + \
                rn*self.opt['RelStopTol']
            edua = np.sqrt(self.Nx)*self.opt['AbsStopTol'] + \
                sn*self.opt['RelStopTol']
        else:
            if rn == 0.0:
                rn = 1.0
            if sn == 0.0:
                sn = 1.0
            r /= rn
            s /= sn
            epri = np.sqrt(self.Nc)*self.opt['AbsStopTol']/rn + \
                self.opt['RelStopTol']
            edua = np.sqrt(self.Nx)*self.opt['AbsStopTol']/sn + \
                self.opt['RelStopTol']

        return r, s, epri, edua
965336
def _Resample(self, stats, target_size):
    """"""
    t_first = stats[0][0]
    t_last = stats[-1][0]
    interval = (t_last - t_first) / target_size

    result = []

    current_t = t_first
    current_v = 0
    i = 0
    while i < len(stats):
      stat_t = stats[i][0]
      stat_v = stats[i][1]
      if stat_t <= (current_t + interval):
        # Always add the last value in an interval to the result.
        current_v = stat_v
        i += 1
      else:
        result.append([current_t + interval, current_v])
        current_t += interval

    result.append([current_t + interval, current_v])
    return result
954688
def resample(self,N):
        
        return rand.random(size=N)*(self.maxval - self.minval) + self.minval
279278
def compute_residuals(self):
        """"""

        r = self.rsdl()
        adapt_tol = self.opt['RelStopTol']

        if self.opt['AutoStop', 'Enabled']:
            adapt_tol = self.tau0 / (1. + self.k)

        return r, adapt_tol