Compute relative residual norms from Hessenberg matrix.

    It is assumed that the initial guess is chosen as zero.
431795
def get_residual_norms(H, self_adjoint=False):
    ''''''
    H = H.copy()
    n_, n = H.shape
    y = numpy.eye(n_, 1, dtype=H.dtype)
    resnorms = [1.]
    for i in range(n_-1):
        G = Givens(H[i:i+2, [i]])
        if self_adjoint:
            H[i:i+2, i:i+3] = G.apply(H[i:i+2, i:i+3])
        else:
            H[i:i+2, i:] = G.apply(H[i:i+2, i:])
        y[i:i+2] = G.apply(y[i:i+2])
        resnorms.append(numpy.abs(y[i+1, 0]))
    if n_ == n:
        resnorms.append(0.)
    return numpy.array(resnorms)
729034
def getSparseTensor(numNonzeros, inputSize, outputSize,
                    onlyPositive=False,
                    fixedRange=1.0/24):
  
  # Initialize weights in the typical fashion.
  w = torch.Tensor(outputSize, inputSize, )

  if onlyPositive:
    w.data.uniform_(0, fixedRange)
  else:
    w.data.uniform_(-fixedRange, fixedRange)

  # Zero out weights for sparse weight matrices
  if numNonzeros < inputSize:
    numZeros = inputSize - numNonzeros

    outputIndices = np.arange(outputSize)
    inputIndices = np.array([np.random.permutation(inputSize)[:numZeros]
                             for _ in outputIndices], dtype=np.long)

    # Create tensor indices for all non-zero weights
    zeroIndices = np.empty((outputSize, numZeros, 2), dtype=np.long)
    zeroIndices[:, :, 0] = outputIndices[:, None]
    zeroIndices[:, :, 1] = inputIndices
    zeroIndices = torch.LongTensor(zeroIndices.reshape(-1, 2))

    zeroWts = (zeroIndices[:, 0], zeroIndices[:, 1])
    w.data[zeroWts] = 0.0

  return w
829424
def apply_matrix_norm(m, v):
    """"""
    (a, b, c, d, e, f) = m
    (p, q) = v
    return (a*p+c*q, b*p+d*q)
435553
public void gaussianElimination(int[][] matrix) {
    assert(matrix != null && matrix.length > 0 && matrix[0].length > 0
           && matrix.length < matrix[0].length);
    int height = matrix.length;
    int width = matrix[0].length;
    for (int i = 0; i < height; i++) {
      boolean pivotFound = false;
      // scan the column for a nonzero pivot and swap it to the diagonal
      for (int j = i; j < height; j++) {
        if (matrix[i][j] != 0) {
          int[] tmp = matrix[i];
          matrix[i] = matrix[j];
          matrix[j] = tmp;
          pivotFound = true;
          break;
        }
      }
      if (!pivotFound) {
        continue;
      }
      int pivot = matrix[i][i];
      for (int j = i; j < width; j++) {
        matrix[i][j] = divide(matrix[i][j], pivot);
      }
      for (int j = i + 1; j < height; j++) {
        int lead = matrix[j][i];
        for (int k = i; k < width; k++) {
          matrix[j][k] = add(matrix[j][k], multiply(lead, matrix[i][k]));
        }
      }
    }
    for (int i = height - 1; i >=0; i--) {
      for (int j = 0; j < i; j++) {
        int lead = matrix[j][i];
        for (int k = i; k < width; k++) {
          matrix[j][k] = add(matrix[j][k], multiply(lead, matrix[i][k]));
        }
      }
    }
  }
1271192
def __basis0(self, xi):
        """"""
        return np.where(np.all([self.knot_vector[:-1] <=  xi,
                                xi < self.knot_vector[1:]],axis=0), 1.0, 0.0)
398767
def _get_kriging_matrix(self, n, n_withdrifts):
        """"""

        xy = np.concatenate((self.X_ADJUSTED[:, np.newaxis],
                             self.Y_ADJUSTED[:, np.newaxis]), axis=1)
        d = cdist(xy, xy, 'euclidean')
        if self.UNBIAS:
            a = np.zeros((n_withdrifts+1, n_withdrifts+1))
        else:
            a = np.zeros((n_withdrifts, n_withdrifts))
        a[:n, :n] = - self.variogram_function(self.variogram_model_parameters, d)
        np.fill_diagonal(a, 0.)

        i = n
        if self.regional_linear_drift:
            a[:n, i] = self.X_ADJUSTED
            a[i, :n] = self.X_ADJUSTED
            i += 1
            a[:n, i] = self.Y_ADJUSTED
            a[i, :n] = self.Y_ADJUSTED
            i += 1
        if self.point_log_drift:
            for well_no in range(self.point_log_array.shape[0]):
                log_dist = np.log(np.sqrt((self.X_ADJUSTED - self.point_log_array[well_no, 0])**2 +
                                          (self.Y_ADJUSTED - self.point_log_array[well_no, 1])**2))
                if np.any(np.isinf(log_dist)):
                    log_dist[np.isinf(log_dist)] = -100.0
                a[:n, i] = - self.point_log_array[well_no, 2] * log_dist
                a[i, :n] = - self.point_log_array[well_no, 2] * log_dist
                i += 1
        if self.external_Z_drift:
            a[:n, i] = self.z_scalars
            a[i, :n] = self.z_scalars
            i += 1
        if self.specified_drift:
            for arr in self.specified_drift_data_arrays:
                a[:n, i] = arr
                a[i, :n] = arr
                i += 1
        if self.functional_drift:
            for func in self.functional_drift_terms:
                a[:n, i] = func(self.X_ADJUSTED, self.Y_ADJUSTED)
                a[i, :n] = func(self.X_ADJUSTED, self.Y_ADJUSTED)
                i += 1
        if i != n_withdrifts:
            warnings.warn("Error in creating kriging matrix. Kriging may fail.",
                          RuntimeWarning)
        if self.UNBIAS:
            a[n_withdrifts, :n] = 1.0
            a[:n, n_withdrifts] = 1.0
            a[n:n_withdrifts + 1, n:n_withdrifts + 1] = 0.0

        return a
1025283
private static double getNorm(double[] w)
  {
    double norm = 0;
    for (int i = 0; i < w.length; i++) {
      norm += w[i]*w[i];
    }
    return norm;
  }
965016
public static double calculateScaledResidual(DoubleMatrix2D A, DoubleMatrix2D X, DoubleMatrix2D B){
		double residual = -Double.MAX_VALUE;
		double niX = Algebra.DEFAULT.normInfinity(X);
		double niB = Algebra.DEFAULT.normInfinity(B);
		if(Double.compare(niX, 0.)==0 && Double.compare(niB, 0.)==0){
			return 0;
		}else{
			double num = Algebra.DEFAULT.normInfinity(Algebra.DEFAULT.mult(A, X).assign(B,	Functions.minus));
			double den = Algebra.DEFAULT.normInfinity(A) * niX + niB;
			residual =  num / den;
			//log.debug("scaled residual: " + residual);
			
			return residual;
		}
	}
1123446
def _init_norm(self, weights):
    """"""
    with tf.variable_scope("init_norm"):
      flat = tf.reshape(weights, [-1, self.layer_depth])
      return tf.reshape(tf.norm(flat, axis=0), (self.layer_depth,))
347171
def _get_noise(self, residuals: np.ndarray) -> float:
        
        return np.mean(np.abs(residuals))
1028362
def _generate_relative_positions_matrix(length_q, length_k,
                                        max_relative_position,
                                        cache=False):
  """"""
  if not cache:
    if length_q == length_k:
      range_vec_q = range_vec_k = tf.range(length_q)
    else:
      range_vec_k = tf.range(length_k)
      range_vec_q = range_vec_k[-length_q:]
    distance_mat = range_vec_k[None, :] - range_vec_q[:, None]
  else:
    distance_mat = tf.expand_dims(tf.range(-length_k+1, 1, 1), 0)
  distance_mat_clipped = tf.clip_by_value(distance_mat, -max_relative_position,
                                          max_relative_position)
  # Shift values to be >= 0. Each integer still uniquely identifies a relative
  # position difference.
  final_mat = distance_mat_clipped + max_relative_position
  return final_mat
497121
def _spectral_norm(self):
        
        w = self.params.get('weight').data(self.ctx)
        w_mat = nd.reshape(w, [w.shape[0], -1])

        _u = self.u.data(self.ctx)
        _v = None

        for _ in range(POWER_ITERATION):
            _v = nd.L2Normalization(nd.dot(_u, w_mat))
            _u = nd.L2Normalization(nd.dot(_v, w_mat.T))

        sigma = nd.sum(nd.dot(_u, w_mat) * _v)
        if sigma == 0.:
            sigma = EPSILON

        with autograd.pause():
            self.u.set_data(_u)

        return w / sigma
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
645953
def relative_idx(colptr, rowidx, snptr, snpar):
    
    
    relptr = matrix(0, (len(snptr),1))
    relidx = matrix(-1, (colptr[-1],1))

    def lfind(a,b):
        i = 0
        ret = +a
        for k in range(len(a)):
            while a[k] != b[i]: i += 1
            ret[k] = i
            i += 1
        return ret
    
    for k in range(len(snpar)):
        p = snpar[k]
        relptr[k+1] = relptr[k]
        if p != -1:
            nk = snptr[k+1] - snptr[k]
            relptr[k+1] += colptr[k+1] - colptr[k] - nk
            relidx[relptr[k]:relptr[k+1]] = lfind(rowidx[colptr[k]+nk:colptr[k+1]], rowidx[colptr[p]:colptr[p+1]])

    return relptr, relidx[:relptr[k+1]]
827941
public static LinearSolverDense<DMatrixRMaj> general(int numRows , int numCols ) {
        if( numRows == numCols )
            return linear(numRows);
        else
            return leastSquares(numRows,numCols);
    }
548195
def _init_tmatrix(self):
        

        if self.radius_type == Scatterer.RADIUS_MAXIMUM:
            # Maximum radius is not directly supported in the original
            # so we convert it to equal volume radius
            radius_type = Scatterer.RADIUS_EQUAL_VOLUME
            radius = self.equal_volume_from_maximum()
        else:
            radius_type = self.radius_type
            radius = self.radius

        self.nmax = pytmatrix.calctmat(radius, radius_type,
            self.wavelength, self.m.real, self.m.imag, self.axis_ratio,
            self.shape, self.ddelt, self.ndgs)
        self._tm_signature = (self.radius, self.radius_type, self.wavelength,
            self.m, self.axis_ratio, self.shape, self.ddelt, self.ndgs)
537614
def normalize_residuals(self, points):
        """"""
        residuals = self.evaluate_residual(points)
        solutions = self.evaluate_solution(points)
        return [resid / soln for resid, soln in zip(residuals, solutions)]