Mix-in class that returns matrix rows for leaky wall condition.
        Qnormal = resfac * (headin - headout)
        Returns matrix part (nunknowns,neq)
        Returns rhs part nunknowns
199049
def equation(self):
        
        mat = np.empty((self.nunknowns, self.model.neq))
        rhs = np.zeros(self.nunknowns)  # Needs to be initialized to zero
        for icp in range(self.ncp):
            istart = icp * self.nlayers
            ieq = 0
            for e in self.model.elementlist:
                if e.nunknowns > 0:
                    qx, qy = e.disvecinflayers(self.xc[icp], self.yc[icp], self.layers)
                    mat[istart:istart + self.nlayers, ieq:ieq + e.nunknowns] = \
                        qx * self.cosnorm[icp] + qy * self.sinnorm[icp] - self.resfac[:, np.newaxis] * \
                        (e.potinflayers(self.xcin[icp], self.ycin[icp], self.layers, aq=self.aq) / self.aq.Tcol[
                            self.layers] - \
                         e.potinflayers(self.xcout[icp], self.ycout[icp], self.layers, aq=self.aq) / self.aq.Tcol[
                             self.layers])
                    ieq += e.nunknowns
                else:
                    qx, qy = e.disveclayers(self.xc[icp], self.yc[icp], self.layers)
                    rhs[istart:istart + self.nlayers] -= qx * self.cosnorm[icp] + qy * self.sinnorm[icp] + self.resfac * \
                                                          (e.potentiallayers(self.xcin[icp], self.ycin[icp], self.layers,
                                                                             aq=self.aq) / self.aq.T[self.layers] -
                                                           e.potentiallayers(self.xcout[icp], self.ycout[icp],
                                                                             self.layers, aq=self.aq) / self.aq.T[
                                                              self.layers])
        return mat, rhs
417706
def lowPassFilter(self, threshold):
        '''
        
        '''
        if not threshold:
            return
        rows, cols = self.img.shape
        tx = int(cols * threshold * 0.25)
        ty = int(rows * threshold * 0.25)
        # upper side
        self.fshift[rows - tx:rows, :] = 0
        # lower side
        self.fshift[0:tx, :] = 0
        # left side
        self.fshift[:, 0:ty] = 0
        # right side
        self.fshift[:, cols - ty:cols] = 0
630546
public static int cudnnBatchNormalizationForwardTraining(
        cudnnHandle handle, 
        int mode, 
        Pointer alpha, /** alpha[0] = result blend factor */
        Pointer beta, /** beta[0] = dest layer blend factor */
        cudnnTensorDescriptor xDesc, 
        Pointer x, /** NxCxHxW */
        cudnnTensorDescriptor yDesc, 
        Pointer y, /** NxCxHxW */
        /**
         * <pre>
         * Shared desc for the next 6 tensors in the argument list.
               Data type to be set as follows:
               type = (typeOf(x) == double) ? double : float
               Dimensions for this descriptor depend on normalization mode
               - Spatial Normalization : tensors are expected to have dims 1xCx1x1
                (normalization is performed across NxHxW)
               - Per-Activation Normalization : tensors are expected to have dims of 1xCxHxW
         * (normalization is performed across N)
         * </pre>
         */
        cudnnTensorDescriptor bnScaleBiasMeanVarDesc, 
        /** 'Gamma' and 'Beta' respectively in Ioffe and Szegedy's paper's notation */
        Pointer bnScale, 
        Pointer bnBias, 
        /**
         * <pre>
         * MUST use factor=1 in the very first call of a complete training cycle.
               Use a factor=1/(1+n) at N-th call to the function to get
               Cumulative Moving Average (CMA) behavior
               CMA[n] = (x[1]+...+x[n])/n
               Since CMA[n+1] = (n*CMA[n]+x[n+1])/(n+1) =
               ((n+1)*CMA[n]-CMA[n])/(n+1) + x[n+1]/(n+1) =
         * CMA[n]*(1-1/(n+1)) + x[n+1]*1/(n+1)
         * </pre>
         */
        double exponentialAverageFactor, 
        /** Used in Training phase only.
               runningMean = newMean*factor + runningMean*(1-factor) */
        Pointer resultRunningMean, 
        /** Output in training mode, input in inference. Is the moving average
               of  variance[x] (factor is applied in the same way as for runningMean) */
        Pointer resultRunningVariance, 
        /** Has to be >= CUDNN_BN_MIN_EPSILON. Should be the same in forward and backward functions. */
        double epsilon, 
        /** Optionally save intermediate results from the forward pass here
               - can be reused to speed up backward pass. NULL if unused */
        Pointer resultSaveMean, 
        Pointer resultSaveInvVariance)
    {
        return checkResult(cudnnBatchNormalizationForwardTrainingNative(handle, mode, alpha, beta, xDesc, x, yDesc, y, bnScaleBiasMeanVarDesc, bnScale, bnBias, exponentialAverageFactor, resultRunningMean, resultRunningVariance, epsilon, resultSaveMean, resultSaveInvVariance));
    }
231345
public static DMatrixRMaj checkZerosLT(DMatrixRMaj A , int numRows , int numCols) {
        if( A == null ) {
            return new DMatrixRMaj(numRows,numCols);
        }  else if( numRows != A.numRows || numCols != A.numCols ) {
            A.reshape(numRows,numCols);
            A.zero();
        } else {
            for( int i = 0; i < A.numRows; i++ ) {
                int index = i*A.numCols;
                int end = index + Math.min(i,A.numCols);;
                while( index < end ) {
                    A.data[index++] = 0;
                }
            }
        }
        return A;
    }
1145937
@Override
    public DMatrixRMaj getR(DMatrixRMaj R, boolean compact) {
        if( R == null ) {
            if( compact ) {
                R = new DMatrixRMaj(minLength,numCols);
            } else
                R = new DMatrixRMaj(numRows,numCols);
        } else {
            if( compact ) {
                R.reshape(minLength,numCols);
            } else {
                R.reshape(numRows,numCols);
            }

            for( int i = 0; i < R.numRows; i++ ) {
                int min = Math.min(i,R.numCols);
                for( int j = 0; j < min; j++ ) {
                    R.set(i,j,0);
                }
            }
        }

        for( int i = 0; i < minLength; i++ ) {
            for( int j = i; j < numCols; j++ ) {
                double val = QR.get(i,j);
                R.set(i,j,val);
            }
        }

        return R;
    }
229214
@Override
    public ZMatrixRMaj getR(ZMatrixRMaj R, boolean compact) {
        if( compact )
            R = UtilDecompositons_ZDRM.checkZerosLT(R,minLength,numCols);
        else
            R = UtilDecompositons_ZDRM.checkZerosLT(R,numRows,numCols);

        for( int i = 0; i < R.numRows; i++ ) {
            for( int j = i; j < R.numCols; j++ ) {
                int index = QR.getIndex(j,i);
                R.set(i,j,QR.data[index],QR.data[index+1]);
            }
        }

        return R;
    }
1144392
@Override
    public DMatrixRMaj getA() {
        if( A.data.length < numRows*numCols ) {
            A = new DMatrixRMaj(numRows,numCols);
        }
        A.reshape(numRows,numCols, false);
        CommonOps_DDRM.mult(Q,R,A);

        return A;
    }
1209275
function getMatrix(data) {
    var matrix = init(data.version);
    fillFinders(matrix);
    fillAlignAndTiming(matrix);
    fillStub(matrix);

    var penalty = Infinity;
    var bestMask = 0;
    for (var mask = 0; mask < 8; mask++) {
        fillData(matrix, data, mask);
        fillReserved(matrix, data.ec_level, mask);
        var p = calculatePenalty(matrix);
        if (p < penalty) {
            penalty = p;
            bestMask = mask;
        }
    }

    fillData(matrix, data, bestMask);
    fillReserved(matrix, data.ec_level, bestMask);

    return matrix.map(function(row) {
        return row.map(function(cell) {
            return cell & 1;
        });
    });
}
205648
def MatrixSolve(a, rhs, adj):
    
    return np.linalg.solve(a if not adj else _adjoint(a), rhs),
888326
public static RQ factorize(Matrix A) {
        return new RQ(A.numRows(), A.numColumns()).factor(new DenseMatrix(A));
    }
588293
public DenseMatrix getQ() {
        int m = qr.nrows();
        int n = qr.ncols();
        DenseMatrix Q = Matrix.zeros(m, n);
        for (int k = n - 1; k >= 0; k--) {
            Q.set(k, k, 1.0);
            for (int j = k; j < n; j++) {
                if (qr.get(k, k) != 0) {
                    double s = 0.0;
                    for (int i = k; i < m; i++) {
                        s += qr.get(i, k) * Q.get(i, j);
                    }
                    s = -s / qr.get(k, k);
                    for (int i = k; i < m; i++) {
                        Q.add(i, j, s * qr.get(i, k));
                    }
                }
            }
        }
        return Q;
    }
1091580
def _transform_coefficients(self, NN, HHw, CCw, ffparm, polycf,
                               any_pwl, npol, nw):
        
        nnw = any_pwl + npol + nw
        M = csr_matrix((ffparm[:, 3], (range(nnw), range(nnw))))
        MR = M * ffparm[:, 2] # FIXME: Possibly column 1.
        HMR = HHw * MR
        MN = M * NN
        HH = MN.T * HHw * MN
        CC = MN.T * (CCw - HMR)
        # Constant term of cost.
        C0 = 1./2. * MR.T * HMR + sum(polycf[:, 2])

        return HH, CC, C0[0]
46134
@Override
    public DMatrixRMaj getQ(DMatrixRMaj Q , boolean compact ) {
        if( compact ) {
            Q = UtilDecompositons_DDRM.checkIdentity(Q,numRows,minLength);
        } else {
            Q = UtilDecompositons_DDRM.checkIdentity(Q,numRows,numRows);
        }

        for( int j = minLength-1; j >= 0; j-- ) {
            double u[] = dataQR[j];

            double vv = u[j];
            u[j] = 1;
            QrHelperFunctions_DDRM.rank1UpdateMultR(Q, u, gammas[j], j, j, numRows, v);
            u[j] = vv;
        }

        return Q;
    }
875481
def mtakfac(m):
  
  u, s, v = msvd(np.asarray(m, dtype=complex))
  f = np.sqrt(u.conj().T @ v.conj())
  w = v @ f
  return w, s
903738
def rmatrixquaternion(q):
    
    assert np.allclose(math.sqrt(np.dot(q,q)), 1.0)

    x, y, z, w = q

    xx = x*x
    xy = x*y
    xz = x*z
    xw = x*w
    yy = y*y
    yz = y*z
    yw = y*w
    zz = z*z
    zw = z*w

    r00 = 1.0 - 2.0 * (yy + zz)
    r01 =       2.0 * (xy - zw)
    r02 =       2.0 * (xz + yw)

    r10 =       2.0 * (xy + zw)
    r11 = 1.0 - 2.0 * (xx + zz)
    r12 =       2.0 * (yz - xw)

    r20 =       2.0 * (xz - yw)
    r21 =       2.0 * (yz + xw)
    r22 = 1.0 - 2.0 * (xx + yy)

    R = np.array([[r00, r01, r02],
               [r10, r11, r12],
               [r20, r21, r22]], float)

    assert np.allclose(np.linalg.det(R), 1.0)
    return R
782016
def dchisq(psr,formbats=False,renormalize=True):
    """"""
    
    if formbats:
        psr.formbats()

    res, err = psr.residuals(removemean=False)[psr.deleted == 0], psr.toaerrs[psr.deleted == 0]

    res -= numpy.sum(res/err**2) / numpy.sum(1/err**2)
    
    # bats already updated by residuals(); skip constant-phase column
    M = psr.designmatrix(updatebats=False,fixunits=True,fixsigns=True)[psr.deleted==0,1:]

    # renormalize design-matrix columns
    if renormalize:
        norm = numpy.sqrt(numpy.sum(M**2,axis=0))
        M /= norm
    else:
        norm = 1.0

    # compute chisq derivative, de-renormalize
    dr = -2 * numpy.dot(M.T,res / (1e-12 * err**2)) * norm
    
    return dr
762913
@Override
    public DMatrixRMaj getQ(DMatrixRMaj Q , boolean compact ) {
        if( compact ) {
            if( Q == null ) {
                Q = CommonOps_DDRM.identity(numRows,minLength);
            } else {
                if( Q.numRows != numRows || Q.numCols != minLength ) {
                    throw new IllegalArgumentException("Unexpected matrix dimension.");
                } else {
                    CommonOps_DDRM.setIdentity(Q);
                }
            }
        } else {
            if( Q == null ) {
                Q = CommonOps_DDRM.identity(numRows);
            } else {
                if( Q.numRows != numRows || Q.numCols != numRows ) {
                    throw new IllegalArgumentException("Unexpected matrix dimension.");
                } else {
                    CommonOps_DDRM.setIdentity(Q);
                }
            }
        }

        for( int j = rank-1; j >= 0; j-- ) {
            double u[] = dataQR[j];

            double vv = u[j];
            u[j] = 1;
            QrHelperFunctions_DDRM.rank1UpdateMultR(Q, u, gammas[j], j, j, numRows, v);
            u[j] = vv;
        }

        return Q;
    }