Python scipy.randn方法代码示例

# 需要导入模块: import scipy [as 别名]
# 或者: from scipy import randn [as 别名]
def haar_measure(n):
    """A Random matrix distributed with the Haar measure.

    For more details, see :cite:`mezzadri2006`.

    Args:
        n (int): matrix size
    Returns:
        array: an nxn random matrix
    """
    z = (sp.randn(n, n) + 1j * sp.randn(n, n)) / np.sqrt(2.0)
    q, r = qr(z)
    d = sp.diagonal(r)
    ph = d / np.abs(d)
    q = np.multiply(q, ph, q)
    return q 

# 需要导入模块: import scipy [as 别名]
# 或者: from scipy import randn [as 别名]
def generate_data(N, S, L):

        # generate genetics
        G = 1.0 * (sp.rand(N, S) < 0.2)
        G -= G.mean(0)
        G /= G.std(0) * sp.sqrt(G.shape[1])

        # generate latent phenotypes
        Zg = sp.dot(G, sp.randn(G.shape[1], L))
        Zn = sp.randn(N, L)

        # generate variance exapleind
        vg = sp.linspace(0.8, 0, L)

        # rescale and sum
        Zg *= sp.sqrt(vg / Zg.var(0))
        Zn *= sp.sqrt((1 - vg) / Zn.var(0))
        Z = Zg + Zn

        return Z, G 

# 需要导入模块: import scipy [as 别名]
# 或者: from scipy import randn [as 别名]
def test_warm_start(self):

        # Big problem
        sp.random.seed(2)
        self.n = 100
        self.m = 200
        self.A = sparse.random(self.m, self.n, density=0.9, format='csc')
        self.l = -sp.rand(self.m) * 2.
        self.u = sp.rand(self.m) * 2.

        P = sparse.random(self.n, self.n, density=0.9)
        self.P = sparse.triu(P.dot(P.T), format='csc')
        self.q = sp.randn(self.n)

        # Setup solver
        self.model = osqp.OSQP()
        self.model.setup(P=self.P, q=self.q, A=self.A, l=self.l, u=self.u,
                         **self.opts)

        # Solve problem with OSQP
        res = self.model.solve()

        # Store optimal values
        x_opt = res.x
        y_opt = res.y
        tot_iter = res.info.iter

        # Warm start with zeros and check if number of iterations is the same
        self.model.warm_start(x=np.zeros(self.n), y=np.zeros(self.m))
        res = self.model.solve()
        self.assertEqual(res.info.iter, tot_iter)

        # Warm start with optimal values and check that number of iter < 10
        self.model.warm_start(x=x_opt, y=y_opt)
        res = self.model.solve()
        self.assertLess(res.info.iter, 10) 

# 需要导入模块: import scipy [as 别名]
# 或者: from scipy import randn [as 别名]
def test_polish_unconstrained(self):

        # Unconstrained QP problem
        sp.random.seed(4)

        self.n = 30
        self.m = 0
        P = sparse.diags(np.random.rand(self.n)) + 0.2*sparse.eye(self.n)
        self.P = P.tocsc()
        self.q = np.random.randn(self.n)
        self.A = sparse.csc_matrix((self.m, self.n))
        self.l = np.array([])
        self.u = np.array([])
        self.model = osqp.OSQP()
        self.model.setup(P=self.P, q=self.q, A=self.A, l=self.l, u=self.u,
                         **self.opts)

        # Solve problem
        res = self.model.solve()

        # Assert close
        nptest.assert_array_almost_equal(
            res.x, np.array([
                -0.61981415, -0.06174194, 0.83824061, -0.0595013, -0.17810828,
                2.90550031, -1.8901713, -1.91191741, -3.73603446, 1.7530356,
                -1.67018181, 3.42221944, 0.61263403, -0.45838347, -0.13194248,
                2.95744794, 5.2902277, -1.42836238, -8.55123842, -0.79093815,
                0.43418189, -0.69323554, 1.15967924, -0.47821898, 3.6108927,
                0.03404309, 0.16322926, -2.17974795, 0.32458796, -1.97553574]))
        nptest.assert_array_almost_equal(res.y, np.array([]))
        nptest.assert_array_almost_equal(res.info.obj_val, -35.020288603855825) 

# 需要导入模块: import scipy [as 别名]
# 或者: from scipy import randn [as 别名]
def test_polish_random(self):

        # Random QP problem
        sp.random.seed(6)

        self.n = 30
        self.m = 50
        Pt = sp.randn(self.n, self.n)
        self.P = sparse.triu(np.dot(Pt.T, Pt), format='csc')
        self.q = sp.randn(self.n)
        self.A = sparse.csc_matrix(sp.randn(self.m, self.n))
        self.l = -3 + sp.randn(self.m)
        self.u = 3 + sp.randn(self.m)
        self.model = osqp.OSQP()
        self.model.setup(P=self.P, q=self.q, A=self.A, l=self.l, u=self.u,
                         **self.opts)

        # Solve problem
        res = self.model.solve()

        # Assert close
        nptest.assert_array_almost_equal(
            res.x, np.array([
                -0.58549607, 0.0030388, -0.07154039, -0.0406463, -0.13349925,
                -0.1354755, -0.17417362, 0.0165324, -0.12213118, -0.10477034,
                -0.51748662, -0.05310921, 0.07862616, 0.53663003, -0.01459859,
                0.40678716, -0.03496123, 0.25722838, 0.06335071, 0.29908295,
                -0.6223218, -0.07614658, -0.3892153, -0.18111635, 0.56301768,
                0.10429917, 0.09821862, -0.30881928, 0.24430531, 0.06597486]))
        nptest.assert_array_almost_equal(
            res.y, np.array([
                0., -2.11407101e-01, 0., 0., 0., 0., 0., 0., 0.,
                0., -3.78854588e-02, 0., -1.58346998e-02, 0., 0.,
                -6.88711599e-02, 0., 0., 0., 0., 0., 0., 0., 0.,
                6.04385132e-01, 0., 0., 0., 0., 0., 0., 0., 0.,
                0., 1.37995470e-01, 0., 0., 0.,  -2.04427802e-02,
                0., -1.32983915e-01, 0., 2.94425952e-02, 0., 0.,
                0., 0., 0., -6.53409219e-02, 0.]))
        nptest.assert_array_almost_equal(res.info.obj_val, -3.262280663471232) 

# 需要导入模块: import scipy [as 别名]
# 或者: from scipy import randn [as 别名]
def test_primal_infeasible_problem(self):

        # Simple QP problem
        sp.random.seed(4)

        self.n = 50
        self.m = 500
        # Generate random Matrices
        Pt = sparse.random(self.n, self.n)
        self.P = sparse.triu(Pt.T.dot(Pt), format='csc')
        self.q = sp.randn(self.n)
        self.A = sparse.random(self.m, self.n).tolil()  # Lil for efficiency
        self.u = 3 + sp.randn(self.m)
        self.l = -3 + sp.randn(self.m)

        # Make random problem primal infeasible
        self.A[int(self.n/2), :] = self.A[int(self.n/2)+1, :]
        self.l[int(self.n/2)] = self.u[int(self.n/2)+1] + 10 * sp.rand()
        self.u[int(self.n/2)] = self.l[int(self.n/2)] + 0.5

        # Convert A to csc
        self.A = self.A.tocsc()

        self.model = osqp.OSQP()
        self.model.setup(P=self.P, q=self.q, A=self.A, l=self.l, u=self.u,
                         **self.opts)

        # Solve problem with OSQP
        res = self.model.solve()

        # Assert close
        self.assertEqual(res.info.status_val,
                         constant('OSQP_PRIMAL_INFEASIBLE')) 

# 需要导入模块: import scipy [as 别名]
# 或者: from scipy import randn [as 别名]
def test(dim=2, clen=10):

    for i in range(clen):
        print("-----------------------------------")
        if dim == 2:
            sz = np.arange(256, 513)
        elif dim == 3:
            sz = np.arange(64, 128)

        np.random.shuffle(sz)
        iscplx = [True, False]
        np.random.shuffle(iscplx)
        if iscplx[0]:
            print("Complex input")
            f = np.array(sc.randn(*sz[:dim])+sc.randn(*sz[:dim])*1j)
        else:
            print("Real input")
            f = np.array(sc.randn(*sz[:dim]))

        isac = [True, False]
        np.random.shuffle(isac)
        if isac[0]:
            print("All curvelets")
        else:
            print("Wavelets at finest scale")

        print(f.shape)

        if dim == 2:
            A = ct.fdct2(f.shape, 6, 32, isac[0], cpx=iscplx[0])
        elif dim == 3:
            A = ct.fdct3(f.shape, 4, 8, isac[0], cpx=iscplx[0])

        x = A.fwd(f)

        if np.allclose(norm(f.flatten(), ord=2), norm(x, ord=2)):
            print('Energy check ok!')
        else:
            print('Problem w energy test')

        fr = A.inv(x)
        if np.allclose(f.flatten(), fr.flatten()):
            print('Inverse check ok!')
        else:
            print('Problem w inverse test')

        print("||f|| = ", norm(f.flatten(), ord=2), f.dtype)
        print("||x|| = ", norm(x, ord=2), x.dtype)
        print("||fr|| = ", norm(fr.flatten(), ord=2), fr.dtype)