Function to compute and plot residuals in units of the uncertainty
182376
def _plot_residuals_to_ax(
    data_all, model_ML, ax, e_unit=u.eV, sed=True, errorbar_opts={}
):
    """"""
    if "group" not in data_all.keys():
        data_all["group"] = np.zeros(len(data_all))

    groups = np.unique(data_all["group"])

    MLf_unit, MLsedf = sed_conversion(model_ML[0], model_ML[1].unit, sed)
    MLene = model_ML[0].to(e_unit)
    MLflux = (model_ML[1] * MLsedf).to(MLf_unit)

    ax.axhline(0, color="k", lw=1, ls="--")

    interp = False
    if data_all["energy"].size != MLene.size or not np.allclose(
        data_all["energy"].value, MLene.value
    ):
        interp = True
        from scipy.interpolate import interp1d

        modelfunc = interp1d(MLene.value, MLflux.value, bounds_error=False)

    for g in groups:
        groupidx = np.where(data_all["group"] == g)

        data = data_all[groupidx]

        notul = ~data["ul"]
        df_unit, dsedf = sed_conversion(data["energy"], data["flux"].unit, sed)
        ene = data["energy"].to(e_unit)
        xerr = u.Quantity((data["energy_error_lo"], data["energy_error_hi"]))
        flux = (data["flux"] * dsedf).to(df_unit)
        dflux = (data["flux_error_lo"] + data["flux_error_hi"]) / 2.0
        dflux = (dflux * dsedf).to(df_unit)[notul]

        if interp:
            difference = flux[notul] - modelfunc(ene[notul]) * flux.unit
        else:
            difference = flux[notul] - MLflux[groupidx][notul]

        # wrap around color and marker cycles
        color = color_cycle[int(g) % len(color_cycle)]
        marker = marker_cycle[int(g) % len(marker_cycle)]

        opts = dict(
            zorder=100,
            marker=marker,
            ls="",
            elinewidth=2,
            capsize=0,
            mec=color,
            mew=0.1,
            ms=6,
            color=color,
        )
        opts.update(errorbar_opts)

        ax.errorbar(
            ene[notul].value,
            (difference / dflux).decompose().value,
            yerr=(dflux / dflux).decompose().value,
            xerr=xerr[:, notul].to(e_unit).value,
            **opts
        )

    from matplotlib.ticker import MaxNLocator

    ax.yaxis.set_major_locator(
        MaxNLocator(5, integer="True", prune="upper", symmetric=True)
    )

    ax.set_ylabel(r"$\Delta\sigma$")
    ax.set_xscale("log")
578443
def error_perturbation_var(C, S):
    

    """
    K = S.shape[0]
    cov = tmatrix_cov(C)
    for i in range(K):
        R = S[i, :, :]
        X[i] = (R[:, :, np.newaxis] * cov * R[:, np.newaxis, :]).sum()
916630
def _estimate_eval_intervals(ritz, indices, indices_remaining,
                                 eps_min=0,
                                 eps_max=0,
                                 eps_res=None):
        '''
        '''
        if len(indices) == 0:
            return utils.Intervals(
                [utils.Interval(mu-resnorm, mu+resnorm)
                 for mu, resnorm in zip(ritz.values, ritz.resnorms)])
        if len(ritz.values) == len(indices):
            raise utils.AssumptionError(
                'selection of all Ritz pairs does not allow estimation.')
        if eps_res is None:
            eps_res = numpy.max(numpy.abs([eps_min, eps_max]))

        # compute quantities for bound
        delta_sel = numpy.linalg.norm(ritz.resnorms[indices], 2)
        delta_non_sel = numpy.linalg.norm(ritz.resnorms[indices_remaining], 2)
        delta = utils.gap(ritz.values[indices],
                          ritz.values[indices_remaining])
        mu_ints = utils.Intervals(
            [utils.Interval(mu+eps_min, mu+eps_max)
             for mu in ritz.values[indices]])
        mu_min = mu_ints.min_abs()

        # check gap assumption
        #if delta_sel + delta_non_sel + eps_max - eps_min >= delta:
        if delta_sel + eps_max - eps_min >= delta:
            raise utils.AssumptionError(
                'delta_sel + delta_non_sel + eps_max - eps_min >= delta'
                + '({0} >= {1}'.format(
                    delta_sel + delta_non_sel + eps_max - eps_min,
                    delta)
                )

        # check assumption on mu_min
        if mu_min == 0:
            raise utils.AssumptionError('mu_min == 0 not allowed')

        # compute eta
        #eta = (delta_sel+eps_res)**2 * (
        #    1/(delta-delta_non_sel-eps_max+eps_min)
        #    + 1/mu_min
        #    )
        #left = - delta_non_sel + eps_min - eta
        #right = delta_non_sel + eps_max + eta
        eta = (delta_sel+eps_res)**2 * (
            1/(delta-eps_max+eps_min)
            + 1/mu_min
            )
        left = eps_min - eta
        right = eps_max + eta

        # return bound
        return utils.Intervals(
            [utils.Interval(mu+left, mu+right)
             for mu in ritz.values[indices_remaining]])
554548
def _build_likelihood(self):
        

        """
        K = self.kern.K(self.X)
        L = tf.cholesky(
            K + tf.eye(tf.shape(self.X)[0], dtype=settings.float_type) * settings.numerics.jitter_level)
        F = tf.matmul(L, self.V) + self.mean_function(self.X)

        return tf.reduce_sum(self.likelihood.logp(F, self.Y))
1081197
def compute_mu_star_confidence(ee, num_trajectories, num_resamples,
                               conf_level):
    '''
    
    '''
    ee_resampled = np.zeros([num_trajectories])
    mu_star_resampled = np.zeros([num_resamples])

    if not 0 < conf_level < 1:
        raise ValueError("Confidence level must be between 0-1.")

    resample_index = np.random.randint(
        len(ee), size=(num_resamples, num_trajectories))
    ee_resampled = ee[resample_index]
    # Compute average of the absolute values over each of the resamples
    mu_star_resampled = np.average(np.abs(ee_resampled), axis=1)

    return norm.ppf(0.5 + conf_level / 2) * mu_star_resampled.std(ddof=1)
1291382
def _compute_magnitude_scaling(self, rup, C):
        
        U, SS, NS, RS = self._get_fault_type_dummy_variables(rup)
        if rup.mag <= C['Mh']:
            return C['e1'] * U + C['e2'] * SS + C['e3'] * NS + C['e4'] * RS + \
                C['e5'] * (rup.mag - C['Mh']) + \
                C['e6'] * (rup.mag - C['Mh']) ** 2
        else:
            return C['e1'] * U + C['e2'] * SS + C['e3'] * NS + C['e4'] * RS + \
                C['e7'] * (rup.mag - C['Mh'])
246129
def residual_norm(A, x, b):
    """"""
    return norm(np.ravel(b) - A*np.ravel(x))
821623
def uncertaintyMap(self, psf, method='convolve', fitParams=None):
        '''
        
        '''

        #ignore background:
        #img[img<0]=0
        ###noise should not influence sharpness uncertainty:
        ##img = median_filter(img, 3)

        # decrease noise in order not to overestimate result:
        img = scaleSignal(self.img, fitParams=fitParams)

        if method == 'convolve':
            #print 'convolve'
            blurred = convolve2d(img, psf, 'same')
            m = abs(img-blurred) / abs(img + blurred)
            m = np.nan_to_num(m)
            m*=self.std**2
            m[m>1]=1
            self.blur_distortion = m
            np.save('blurred', blurred)
            return m
        else:
            restored = unsupervised_wiener(img, psf)[0]
            m = abs(img-restored) / abs(img + restored)
            m = np.nan_to_num(m)
            m*=self.std**2
            m[m>1]=1
            self.blur_distortion = m
            return m, restored
518376
def inference_likelihood(self, kern, X, Z, likelihood, Y):
        

        num_data, output_dim = Y.shape
        input_dim = Z.shape[0]
        if self.mpi_comm != None:
            from mpi4py import MPI
            num_data_all = np.array(num_data,dtype=np.int32)
            self.mpi_comm.Allreduce([np.int32(num_data), MPI.INT], [num_data_all, MPI.INT])
            num_data = num_data_all

        if isinstance(X, VariationalPosterior):
            uncertain_inputs = True
        else:
            uncertain_inputs = False

        #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
        if het_noise:
            self.batchsize = 1

        psi0_full, psi1Y_full, psi2_full, YRY_full = self.gatherPsiStat(kern, X, Z, Y, beta, uncertain_inputs)

        #======================================================================
        # Compute Common Components
        #======================================================================

        Kmm = kern.K(Z).copy()
        diag.add(Kmm, self.const_jitter)
        if not np.isfinite(Kmm).all():
            print(Kmm)
        Lm = jitchol(Kmm)
        LmInv = dtrtri(Lm)

        LmInvPsi2LmInvT = LmInv.dot(psi2_full.dot(LmInv.T))
        Lambda = np.eye(Kmm.shape[0])+LmInvPsi2LmInvT
        LL = jitchol(Lambda)
        LLInv = dtrtri(LL)
        logdet_L = 2.*np.sum(np.log(np.diag(LL)))
        LmLLInv = LLInv.dot(LmInv)
        
        b  = psi1Y_full.dot(LmLLInv.T)
        bbt = np.square(b).sum()
        v = b.dot(LmLLInv).T
        LLinvPsi1TYYTPsi1LLinvT = tdot(b.T)
        
        tmp = -LLInv.T.dot(LLinvPsi1TYYTPsi1LLinvT+output_dim*np.eye(input_dim)).dot(LLInv)
        dL_dpsi2R = LmInv.T.dot(tmp+output_dim*np.eye(input_dim)).dot(LmInv)/2.
        
        # Cache intermediate results
        self.midRes['dL_dpsi2R'] = dL_dpsi2R
        self.midRes['v'] = v

        #======================================================================
        # Compute log-likelihood
        #======================================================================
        if het_noise:
            logL_R = -np.log(beta).sum()
        else:
            logL_R = -num_data*np.log(beta)
        logL = -(output_dim*(num_data*log_2_pi+logL_R+psi0_full-np.trace(LmInvPsi2LmInvT))+YRY_full-bbt)/2.-output_dim*logdet_L/2.

        #======================================================================
        # Compute dL_dKmm
        #======================================================================

        dL_dKmm =  dL_dpsi2R - output_dim*LmInv.T.dot(LmInvPsi2LmInvT).dot(LmInv)/2.

        #======================================================================
        # Compute the Posterior distribution of inducing points p(u|Y)
        #======================================================================

        if not self.Y_speedup or het_noise:
            wd_inv = backsub_both_sides(Lm, np.eye(input_dim)- backsub_both_sides(LL, np.identity(input_dim), transpose='left'), transpose='left')
            post = Posterior(woodbury_inv=wd_inv, woodbury_vector=v, K=Kmm, mean=None, cov=None, K_chol=Lm)
        else:
            post = None

        #======================================================================
        # Compute dL_dthetaL for uncertian input and non-heter noise
        #======================================================================

        if not het_noise:
            dL_dthetaL = (YRY_full*beta + beta*output_dim*psi0_full - num_data*output_dim*beta)/2. - beta*(dL_dpsi2R*psi2_full).sum() - beta*(v.T*psi1Y_full).sum()
            self.midRes['dL_dthetaL'] = dL_dthetaL

        return logL, dL_dKmm, post
675501
private void computeGradientAndHessian(DMatrixRMaj param  )
    {
        // residuals = f(x) - y
        function.compute(param, residuals);

        computeNumericalJacobian(param,jacobian);

        CommonOps_DDRM.multTransA(jacobian, residuals, g);
        CommonOps_DDRM.multTransA(jacobian, jacobian,  H);

        CommonOps_DDRM.extractDiag(H,Hdiag);
    }
307831
def Kdiag(self, X):
        """"""
        Kdiag = np.zeros(X.shape[0])
        ly=1/self.lengthscale_Y
        lu=np.sqrt(3)/self.lengthscale_U

        Vu = self.variance_U
        Vy=self.variance_Y

        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

        slices = index_to_slices(X[:,-1])

        for i, ss1 in enumerate(slices):
            for s1 in ss1:
                if i==0:
                    Kdiag[s1]+= self.variance_U + self.ubias
                elif i==1:
                    Kdiag[s1]+= Vu*Vy*(k1+k2+k3)
                else:
                    raise ValueError("invalid input/output index")
        #Kdiag[slices[0][0]]+= self.variance_U   #matern32 diag
        #Kdiag[slices[1][0]]+= self.variance_U*self.variance_Y*(k1+k2+k3)  #  diag
        return Kdiag
307513
def compute(self):
        """"""
        self._compute_window_size()
        smooth = []
        residual = []

        x, y = self.x, self.y

        # step through x and y data with a window window_size wide.
        self._update_values_in_window()
        self._update_mean_in_window()
        self._update_variance_in_window()
        for i, (xi, yi) in enumerate(zip(x, y)):
            if ((i - self._neighbors_on_each_side) > 0.0 and
                    (i + self._neighbors_on_each_side) < len(x)):
                self._advance_window()
            smooth_here = self._compute_smooth_during_construction(xi)
            residual_here = self._compute_cross_validated_residual_here(xi, yi, smooth_here)
            smooth.append(smooth_here)
            residual.append(residual_here)

        self._store_unsorted_results(smooth, residual)
990102
def root_mean_squared_error(pred:Tensor, targ:Tensor)->Rank0Tensor:
    ""
    pred,targ = flatten_check(pred,targ)
    return torch.sqrt(F.mse_loss(pred, targ))
95046
def _assess_approximation(self, boundary_points, derivs, funcs, nodes, problem):
        
        interior_residuals = self._compute_interior_residuals(derivs, funcs,
                                                              nodes, problem)
        boundary_residuals = self._compute_boundary_residuals(boundary_points,
                                                              funcs, problem)
        return np.hstack(interior_residuals + boundary_residuals)
766984
def _compute_cross_validated_residual_here(self, xi, yi, smooth_here):
        
        denom = (1.0 - 1.0 / self.window_size -
                 (xi - self._mean_x_in_window) ** 2 /
                 self._variance_in_window)
        if denom == 0.0:
            # can happen  with small data sets
            return 1.0
        return abs((yi - smooth_here) / denom)
923909
def err(self):
        
        # Initiallize error to nan
        err = odict(self.mle())
        err.update([(k,np.nan*np.ones(2)) for k in err.keys()])

        # Find the maximum likelihood
        a = self.log_likelihood_sparse_array
        j,k = np.unravel_index(a.argmax(),a.shape)

        self.loglike.set_params(distance_modulus=self.distance_modulus_array[j],
                                lon=self.roi.pixels_target.lon[k],
                                lat=self.roi.pixels_target.lat[k])
        self.loglike.sync_params()

        # Find the error in richness, starting at maximum
        lo,hi = np.array(self.loglike.richness_interval())
        err['richness'] = np.array([lo,hi])

        # ADW: This is a rough estimate of the distance uncertainty 
        # hacked to keep the maximum distance modulus on a grid index

        # This is a hack to get the confidence interval to play nice...
        # Require at least three points.
        if (a[:,k]>0).sum() >= 3:
            parabola = Parabola(np.insert(self.distance_modulus_array,0,0.), 
                                np.insert(a[:,k],0,0.) )
            lo,hi = np.array(parabola.confidenceInterval())
            err['distance_modulus'] = self.distance_modulus_array[j] + (hi-lo)/2.*np.array([-1.,1.])

        # ADW: Could estimate lon and lat from the grid.
        # This is just concept right now...
        if (a[j,:]>0).sum() >= 10:
            delta_ts = 2*(a[j,k] - a[j,:])
            pix = np.where(a[j,:][delta_ts < 2.71])[0]
            lons = self.roi.pixels_target.lon[pix]
            lats = self.roi.pixels_target.lat[pix]
            err['lon'] = np.array([ np.min(lons),np.max(lons)])
            err['lat'] = np.array([ np.min(lats),np.max(lats)])

        return err
675129
private void computeDerivatives(double sigma) {
		this.Lx2 = new float[width][height];
		this.Ly2 = new float[width][height];
		this.Lxy = new float[width][height];

		float[][][] grad = new float[width][height][];
		for (int y = 0; y < height; y++)
			for (int x = 0; x < width; x++)
				grad[x][y] = sobel(x, y);

		int radius = (int) (2 * sigma);
		int window = 1 + 2 * radius;
		float[][] gaussian = new float[window][window];
		for (int j = -radius; j <= radius; j++)
			for (int i = -radius; i <= radius; i++)
				gaussian[i + radius][j + radius] = (float) gaussian(i, j, sigma);

		for (int y = 0; y < height; y++) {
			for (int x = 0; x < width; x++) {

				for (int dy = -radius; dy <= radius; dy++) {
					for (int dx = -radius; dx <= radius; dx++) {
						int xk = x + dx;
						int yk = y + dy;
						if (xk < 0 || xk >= width)
							continue;
						if (yk < 0 || yk >= height)
							continue;

						double gw = gaussian[dx + radius][dy + radius];

						this.Lx2[x][y] += gw * grad[xk][yk][0]
								* grad[xk][yk][0];
						this.Ly2[x][y] += gw * grad[xk][yk][1]
								* grad[xk][yk][1];
						this.Lxy[x][y] += gw * grad[xk][yk][0]
								* grad[xk][yk][1];
					}
				}
			}
		}
	}