Constructs the matrix A and the vector b from a timeseries of toy
    action-angles AA to solve for the vector x = (J_0,J_1,J_2,S...) where
    x contains all Fourier components of the generating function with |n|<N_max
2344
def solver(AA, N_max, symNx = 2, throw_out_modes=False):
    

    # Find all integer component n_vectors which lie within sphere of radius N_max
    # Here we have assumed that the potential is symmetric x->-x, y->-y, z->-z
    # This can be relaxed by changing symN to 1
    # Additionally due to time reversal symmetry S_n = -S_-n so we only consider
    # "half" of the n-vector-space

    angs = unroll_angles(AA.T[3:].T,np.ones(3))

    symNz = 2
    NNx = range(-N_max, N_max+1, symNx)
    NNy = range(-N_max, N_max+1, symNz)
    NNz = range(-N_max, N_max+1, symNz)
    n_vectors = np.array([[i,j,k] for (i,j,k) in product(NNx,NNy,NNz)
                          if(not(i==0 and j==0 and k==0)            # exclude zero vector
                             and (k>0                               # northern hemisphere
                                  or (k==0 and j>0)                 # half of x-y plane
                                  or (k==0 and j==0 and i>0))       # half of x axis
                             and np.sqrt(i*i+j*j+k*k)<=N_max)])     # inside sphere

    xxx = check_each_direction(n_vectors,angs)

    if(throw_out_modes):
        n_vectors = np.delete(n_vectors,check_each_direction(n_vectors,angs),axis=0)

    n = len(n_vectors)+3
    b = np.zeros(shape=(n, ))
    a = np.zeros(shape=(n,n))

    a[:3,:3]=len(AA)*np.identity(3)

    for i in AA:
        a[:3,3:]+=2.*n_vectors.T[:3]*np.cos(np.dot(n_vectors,i[3:]))
        a[3:,3:]+=4.*np.dot(n_vectors,n_vectors.T)*np.outer(np.cos(np.dot(n_vectors,i[3:])),np.cos(np.dot(n_vectors,i[3:])))
        b[:3]+=i[:3]
        b[3:]+=2.*np.dot(n_vectors,i[:3])*np.cos(np.dot(n_vectors,i[3:]))

    a[3:,:3]=a[:3,3:].T

    return np.array(solve(a,b)), n_vectors
1304977
def spin1y_from_xi1_phi_a_phi_s(xi1, phi_a, phi_s):
    
    phi1 = phi1_from_phi_a_phi_s(phi_s, phi_a)
    return xi1 * numpy.sin(phi1)
707982
def F_t(X, Y, S, M_E, E, m0, rho):
    '''  '''
    r = X.shape[1]
    out1 = (((np.dot(np.dot(X, S), Y.T) - M_E) * E)**2).sum() / 2
    out2 = rho * G(Y, m0, r)
    out3 = rho * G(X, m0, r)

    return out1 + out2 + out3
1223802
def rpy_matrix(roll, pitch, yaw):
    """"""
    return np.dot(Rz_matrix(yaw), np.dot(Ry_matrix(pitch), Rx_matrix(roll)))
874405
public static ZMatrixRMaj householderVector(ZMatrixRMaj x ) {
        ZMatrixRMaj u = x.copy();

        double max = CommonOps_ZDRM.elementMaxAbs(u);

        CommonOps_ZDRM.elementDivide(u, max, 0, u);

        double nx = NormOps_ZDRM.normF(u);
        Complex_F64 c = new Complex_F64();
        u.get(0,0,c);

        double realTau,imagTau;

        if( c.getMagnitude() == 0 ) {
            realTau = nx;
            imagTau = 0;
        } else {
            realTau = c.real/c.getMagnitude()*nx;
            imagTau = c.imaginary/c.getMagnitude()*nx;
        }

        u.set(0,0,c.real + realTau,c.imaginary + imagTau);
        CommonOps_ZDRM.elementDivide(u,u.getReal(0,0),u.getImag(0,0),u);

        return u;
    }
480917
public void process( DMatrixRMaj F , List<AssociatedPair> observations ,
						 int width , int height ) {

		int centerX = width/2;
		int centerY = height/2;

		MultiViewOps.extractEpipoles(F, epipole1, epipole2);
		checkEpipoleInside(width, height);


		// compute the transform H which will send epipole2 to infinity
		SimpleMatrix R = rotateEpipole(epipole2,centerX,centerY);
		SimpleMatrix T = translateToOrigin(centerX,centerY);
		SimpleMatrix G = computeG(epipole2,centerX,centerY);

		SimpleMatrix H = G.mult(R).mult(T);

		//Find the two matching transforms
		SimpleMatrix Hzero = computeHZero(F,epipole2,H);

		SimpleMatrix Ha = computeAffineH(observations,H.getDDRM(),Hzero.getDDRM());

		rect1.set(Ha.mult(Hzero).getDDRM());
		rect2.set(H.getDDRM());
	}
1238165
def trt_pmf(matrices):
    
    ntrts, nmags, ndists, nlons, nlats, neps = matrices.shape
    pmf = numpy.zeros(ntrts)
    for t in range(ntrts):
        pmf[t] = 1. - numpy.prod(
            [1. - matrices[t, i, j, k, l, m]
             for i in range(nmags)
             for j in range(ndists)
             for k in range(nlons)
             for l in range(nlats)
             for m in range(neps)])
    return pmf
814682
def xstep(self):
        """

        self.YU[:] = self.Y - self.U
        Zf = sl.rfftn(self.YU, None, self.cri.axisN)
        ZfQ = sl.dot(self.Q.T, Zf, axis=self.cri.axisC)
        b = self.DSfBQ + self.rho * ZfQ

        Xh = sl.solvedbi_sm(self.gDf, self.rho, b, self.c,
                            axis=self.cri.axisM)
        self.Xf[:] = sl.dot(self.Q, Xh, axis=self.cri.axisC)
        self.X = sl.irfftn(self.Xf, self.cri.Nv, self.cri.axisN)

        if self.opt['LinSolveCheck']:
            DDXf = np.conj(self.Df) *  sl.inner(self.Df, self.Xf,
                                                axis=self.cri.axisM)
            DDXfBB = sl.dot(self.B.T.dot(self.B), DDXf, axis=self.cri.axisC)
            ax = DDXfBB + self.rho * self.Xf
            b = sl.dot(self.B.T, self.DSf, axis=self.cri.axisC) + \
                self.rho * Zf
            self.xrrs = sl.rrs(ax, b)
        else:
            self.xrrs = None
41749
def matrix(self, angle):
        """"""
        _rot_mat = {
            'x': self._x_rot,
            'y': self._y_rot,
            'z': self._z_rot
        }
        return _rot_mat[self.axis](angle)
430929
public Matrix4x3f transformAab(float minX, float minY, float minZ, float maxX, float maxY, float maxZ, Vector3f outMin, Vector3f outMax) {
        float xax = m00 * minX, xay = m01 * minX, xaz = m02 * minX;
        float xbx = m00 * maxX, xby = m01 * maxX, xbz = m02 * maxX;
        float yax = m10 * minY, yay = m11 * minY, yaz = m12 * minY;
        float ybx = m10 * maxY, yby = m11 * maxY, ybz = m12 * maxY;
        float zax = m20 * minZ, zay = m21 * minZ, zaz = m22 * minZ;
        float zbx = m20 * maxZ, zby = m21 * maxZ, zbz = m22 * maxZ;
        float xminx, xminy, xminz, yminx, yminy, yminz, zminx, zminy, zminz;
        float xmaxx, xmaxy, xmaxz, ymaxx, ymaxy, ymaxz, zmaxx, zmaxy, zmaxz;
        if (xax < xbx) {
            xminx = xax;
            xmaxx = xbx;
        } else {
            xminx = xbx;
            xmaxx = xax;
        }
        if (xay < xby) {
            xminy = xay;
            xmaxy = xby;
        } else {
            xminy = xby;
            xmaxy = xay;
        }
        if (xaz < xbz) {
            xminz = xaz;
            xmaxz = xbz;
        } else {
            xminz = xbz;
            xmaxz = xaz;
        }
        if (yax < ybx) {
            yminx = yax;
            ymaxx = ybx;
        } else {
            yminx = ybx;
            ymaxx = yax;
        }
        if (yay < yby) {
            yminy = yay;
            ymaxy = yby;
        } else {
            yminy = yby;
            ymaxy = yay;
        }
        if (yaz < ybz) {
            yminz = yaz;
            ymaxz = ybz;
        } else {
            yminz = ybz;
            ymaxz = yaz;
        }
        if (zax < zbx) {
            zminx = zax;
            zmaxx = zbx;
        } else {
            zminx = zbx;
            zmaxx = zax;
        }
        if (zay < zby) {
            zminy = zay;
            zmaxy = zby;
        } else {
            zminy = zby;
            zmaxy = zay;
        }
        if (zaz < zbz) {
            zminz = zaz;
            zmaxz = zbz;
        } else {
            zminz = zbz;
            zmaxz = zaz;
        }
        outMin.x = xminx + yminx + zminx + m30;
        outMin.y = xminy + yminy + zminy + m31;
        outMin.z = xminz + yminz + zminz + m32;
        outMax.x = xmaxx + ymaxx + zmaxx + m30;
        outMax.y = xmaxy + ymaxy + zmaxy + m31;
        outMax.z = xmaxz + ymaxz + zmaxz + m32;
        return this;
    }
565991
def solve(self, lam):
        ''''''
        s = weighted_graphtf(self.nnodes, self.y, self.weights, lam,
                             self.Dk.shape[0], self.Dk.shape[1], self.Dk.nnz,
                             self.Dk.row.astype('int32'), self.Dk.col.astype('int32'), self.Dk.data.astype('double'),
                             self.maxsteps, self.converge,
                             self.beta, self.u)
        self.steps.append(s)
        return self.beta
222120
def fourier_series(x, *a):
    
    output = 0
    output += a[0]/2
    w = a[1]
    for n in range(2, len(a), 2):
        n_ = n/2
        val1 = a[n]
        val2 = a[n+1]
        output += val1*np.sin(n_*x*w)
        output += val2*np.cos(n_*x*w)
    return output
690707
def Ax(self):
        
        #  0         -self.a3   self.a2
        #  self.a3    0        -self.a1
        # -self.a2    self.a1   0        
        m = np.zeros((len(self.a1), 3, 3))
        m[:, 0, 1] = -self.a3
        m[:, 0, 2] = +self.a2
        m[:, 1, 0] = +self.a3
        m[:, 1, 2] = -self.a1
        m[:, 2, 0] = -self.a2
        m[:, 2, 1] = +self.a1        
        return m
740835
def Itax(mt, x, t):
    
    return (Sx(mt, x + 1) - Sx(mt, x + t + 1)) / Dx(mt, x)
470058
def get_matrix(theta, phi, psi, dx, dy, dz,
                   matrix=np.zeros((4, 4), dtype=DTYPE),
                   angles=np.zeros(3, dtype=DTYPE)):
        
        # NB!: matrix and angles by default are being overwritten on each call
        # thus, only created once at compile time.

        angles[0] = theta
        angles[1] = phi
        angles[2] = psi

        cx, cy, cz = np.cos(angles)
        sx, sy, sz = np.sin(angles)

        rotation = matrix[:3, :3]
        rotation.flat = (cx * cz - sx * cy * sz,
                         cx * sz + sx * cy * cz, sx * sy,
                         -sx * cz - cx * cy * sz,
                         -sx * sz + cx * cy * cz, cx * sy,
                         sy * sz,
                         -sy * cz, cy)

        # Translation component
        matrix[:3, 3] = dx, dy, dz
        matrix[3, 3] = 1.
        return matrix
997427
public Vector4f hermite(Vector4fc t0, Vector4fc v1, Vector4fc t1, float t, Vector4f dest) {
        float t2 = t * t;
        float t3 = t2 * t;
        dest.x = (x + x - v1.x() - v1.x() + t1.x() + t0.x()) * t3 + (3.0f * v1.x() - 3.0f * x - t0.x() - t0.x() - t1.x()) * t2 + x * t + x;
        dest.y = (y + y - v1.y() - v1.y() + t1.y() + t0.y()) * t3 + (3.0f * v1.y() - 3.0f * y - t0.y() - t0.y() - t1.y()) * t2 + y * t + y;
        dest.z = (z + z - v1.z() - v1.z() + t1.z() + t0.z()) * t3 + (3.0f * v1.z() - 3.0f * z - t0.z() - t0.z() - t1.z()) * t2 + z * t + z;
        dest.w = (w + w - v1.w() - v1.w() + t1.w() + t0.w()) * t3 + (3.0f * v1.w() - 3.0f * w - t0.w() - t0.w() - t1.w()) * t2 + w * t + w;
        return dest;
    }
905814
def _fftconv(a, b, axes=(0, 1)):
    """"""

    if cp.isrealobj(a) and cp.isrealobj(b):
        fft = cp.fft.rfftn
        ifft = cp.fft.irfftn
    else:
        fft = cp.fft.fftn
        ifft = cp.fft.ifftn
    dims = cp.maximum(cp.asarray([a.shape[i] for i in axes]),
                      cp.asarray([b.shape[i] for i in axes]))
    dims = [int(d) for d in dims]
    af = fft(a, dims, axes)
    bf = fft(b, dims, axes)
    return ifft(af * bf, dims, axes)