1#!/usr/bin/python
2# -*- coding: utf-8 -*-
3
4from mpmath import mp
5from mpmath import libmp
6
7xrange = libmp.backend.xrange
8
9def run_eigsy(A, verbose = False):
10    if verbose:
11        print("original matrix:\n", str(A))
12
13    D, Q = mp.eigsy(A)
14    B = Q * mp.diag(D) * Q.transpose()
15    C = A - B
16    E = Q * Q.transpose() - mp.eye(A.rows)
17
18    if verbose:
19        print("eigenvalues:\n", D)
20        print("eigenvectors:\n", Q)
21
22    NC = mp.mnorm(C)
23    NE = mp.mnorm(E)
24
25    if verbose:
26        print("difference:", NC, "\n", C, "\n")
27        print("difference:", NE, "\n", E, "\n")
28
29    eps = mp.exp( 0.8 * mp.log(mp.eps))
30
31    assert NC < eps
32    assert NE < eps
33
34    return NC
35
36def run_eighe(A, verbose = False):
37    if verbose:
38        print("original matrix:\n", str(A))
39
40    D, Q = mp.eighe(A)
41    B = Q * mp.diag(D) * Q.transpose_conj()
42    C = A - B
43    E = Q * Q.transpose_conj() - mp.eye(A.rows)
44
45    if verbose:
46        print("eigenvalues:\n", D)
47        print("eigenvectors:\n", Q)
48
49    NC = mp.mnorm(C)
50    NE = mp.mnorm(E)
51
52    if verbose:
53        print("difference:", NC, "\n", C, "\n")
54        print("difference:", NE, "\n", E, "\n")
55
56    eps = mp.exp( 0.8 * mp.log(mp.eps))
57
58    assert NC < eps
59    assert NE < eps
60
61    return NC
62
63def run_svd_r(A, full_matrices = False, verbose = True):
64
65    m, n = A.rows, A.cols
66
67    eps = mp.exp(0.8 * mp.log(mp.eps))
68
69    if verbose:
70        print("original matrix:\n", str(A))
71        print("full", full_matrices)
72
73    U, S0, V = mp.svd_r(A, full_matrices = full_matrices)
74
75    S = mp.zeros(U.cols, V.rows)
76    for j in xrange(min(m, n)):
77        S[j,j] = S0[j]
78
79    if verbose:
80        print("U:\n", str(U))
81        print("S:\n", str(S0))
82        print("V:\n", str(V))
83
84    C = U * S * V - A
85    err = mp.mnorm(C)
86    if verbose:
87        print("C\n", str(C), "\n", err)
88    assert err < eps
89
90    D = V * V.transpose() - mp.eye(V.rows)
91    err = mp.mnorm(D)
92    if verbose:
93        print("D:\n", str(D), "\n", err)
94    assert err < eps
95
96    E = U.transpose() * U - mp.eye(U.cols)
97    err = mp.mnorm(E)
98    if verbose:
99        print("E:\n", str(E), "\n", err)
100    assert err < eps
101
102def run_svd_c(A, full_matrices = False, verbose = True):
103
104    m, n = A.rows, A.cols
105
106    eps = mp.exp(0.8 * mp.log(mp.eps))
107
108    if verbose:
109        print("original matrix:\n", str(A))
110        print("full", full_matrices)
111
112    U, S0, V = mp.svd_c(A, full_matrices = full_matrices)
113
114    S = mp.zeros(U.cols, V.rows)
115    for j in xrange(min(m, n)):
116        S[j,j] = S0[j]
117
118    if verbose:
119        print("U:\n", str(U))
120        print("S:\n", str(S0))
121        print("V:\n", str(V))
122
123    C = U * S * V - A
124    err = mp.mnorm(C)
125    if verbose:
126        print("C\n", str(C), "\n", err)
127    assert err  < eps
128
129    D = V * V.transpose_conj() - mp.eye(V.rows)
130    err = mp.mnorm(D)
131    if verbose:
132        print("D:\n", str(D), "\n", err)
133    assert err < eps
134
135    E = U.transpose_conj() * U - mp.eye(U.cols)
136    err = mp.mnorm(E)
137    if verbose:
138        print("E:\n", str(E), "\n", err)
139    assert err < eps
140
141def run_gauss(qtype, a, b):
142    eps = 1e-5
143
144    d, e = mp.gauss_quadrature(len(a), qtype)
145    d -= mp.matrix(a)
146    e -= mp.matrix(b)
147
148    assert mp.mnorm(d) < eps
149    assert mp.mnorm(e) < eps
150
151def irandmatrix(n, range = 10):
152    """
153    random matrix with integer entries
154    """
155    A = mp.matrix(n, n)
156    for i in xrange(n):
157        for j in xrange(n):
158            A[i,j]=int( (2 * mp.rand() - 1) * range)
159    return A
160
161#######################
162
163def test_eighe_fixed_matrix():
164    A = mp.matrix([[2, 3], [3, 5]])
165    run_eigsy(A)
166    run_eighe(A)
167
168    A = mp.matrix([[7, -11], [-11, 13]])
169    run_eigsy(A)
170    run_eighe(A)
171
172    A = mp.matrix([[2, 11, 7], [11, 3, 13], [7, 13, 5]])
173    run_eigsy(A)
174    run_eighe(A)
175
176    A = mp.matrix([[2, 0, 7], [0, 3, 1], [7, 1, 5]])
177    run_eigsy(A)
178    run_eighe(A)
179
180    #
181
182    A = mp.matrix([[2, 3+7j], [3-7j, 5]])
183    run_eighe(A)
184
185    A = mp.matrix([[2, -11j, 0], [+11j, 3, 29j], [0, -29j, 5]])
186    run_eighe(A)
187
188    A = mp.matrix([[2, 11 + 17j, 7 + 19j], [11 - 17j, 3, -13 + 23j], [7 - 19j, -13 - 23j, 5]])
189    run_eighe(A)
190
191def test_eigsy_randmatrix():
192    N = 5
193
194    for a in xrange(10):
195        A = 2 * mp.randmatrix(N, N) - 1
196
197        for i in xrange(0, N):
198            for j in xrange(i + 1, N):
199                A[j,i] = A[i,j]
200
201        run_eigsy(A)
202
203def test_eighe_randmatrix():
204    N = 5
205
206    for a in xrange(10):
207        A = (2 * mp.randmatrix(N, N) - 1) + 1j * (2 * mp.randmatrix(N, N) - 1)
208
209        for i in xrange(0, N):
210            A[i,i] = mp.re(A[i,i])
211            for j in xrange(i + 1, N):
212                A[j,i] = mp.conj(A[i,j])
213
214        run_eighe(A)
215
216def test_eigsy_irandmatrix():
217    N = 4
218    R = 4
219
220    for a in xrange(10):
221        A=irandmatrix(N, R)
222
223        for i in xrange(0, N):
224            for j in xrange(i + 1, N):
225                A[j,i] = A[i,j]
226
227        run_eigsy(A)
228
229def test_eighe_irandmatrix():
230    N = 4
231    R = 4
232
233    for a in xrange(10):
234        A=irandmatrix(N, R) + 1j * irandmatrix(N, R)
235
236        for i in xrange(0, N):
237            A[i,i] = mp.re(A[i,i])
238            for j in xrange(i + 1, N):
239                A[j,i] = mp.conj(A[i,j])
240
241        run_eighe(A)
242
243def test_svd_r_rand():
244    for i in xrange(5):
245        full = mp.rand() > 0.5
246        m = 1 + int(mp.rand() * 10)
247        n = 1 + int(mp.rand() * 10)
248        A = 2 * mp.randmatrix(m, n) - 1
249        if mp.rand() > 0.5:
250            A *= 10
251            for x in xrange(m):
252                for y in xrange(n):
253                    A[x,y]=int(A[x,y])
254
255        run_svd_r(A, full_matrices = full, verbose = False)
256
257def test_svd_c_rand():
258    for i in xrange(5):
259        full = mp.rand() > 0.5
260        m = 1 + int(mp.rand() * 10)
261        n = 1 + int(mp.rand() * 10)
262        A = (2 * mp.randmatrix(m, n) - 1) + 1j * (2 * mp.randmatrix(m, n) - 1)
263        if mp.rand() > 0.5:
264            A *= 10
265            for x in xrange(m):
266                for y in xrange(n):
267                    A[x,y]=int(mp.re(A[x,y])) + 1j * int(mp.im(A[x,y]))
268
269        run_svd_c(A, full_matrices=full, verbose=False)
270
271def test_svd_test_case():
272    # a test case from Golub and Reinsch
273    #  (see wilkinson/reinsch: handbook for auto. comp., vol ii-linear algebra, 134-151(1971).)
274
275    eps = mp.exp(0.8 * mp.log(mp.eps))
276
277    a = [[22, 10,  2,   3,  7],
278         [14,  7, 10,   0,  8],
279         [-1, 13, -1, -11,  3],
280         [-3, -2, 13,  -2,  4],
281         [ 9,  8,  1,  -2,  4],
282         [ 9,  1, -7,   5, -1],
283         [ 2, -6,  6,   5,  1],
284         [ 4,  5,  0,  -2,  2]]
285
286    a = mp.matrix(a)
287    b = mp.matrix([mp.sqrt(1248), 20, mp.sqrt(384), 0, 0])
288
289    S = mp.svd_r(a, compute_uv = False)
290    S -= b
291    assert mp.mnorm(S) < eps
292
293    S = mp.svd_c(a, compute_uv = False)
294    S -= b
295    assert mp.mnorm(S) < eps
296
297
298def test_gauss_quadrature_static():
299    a = [-0.57735027,  0.57735027]
300    b = [ 1,  1]
301    run_gauss("legendre", a , b)
302
303    a = [ -0.906179846,  -0.538469310,   0,           0.538469310,   0.906179846]
304    b = [  0.23692689,    0.47862867,    0.56888889,  0.47862867,    0.23692689]
305    run_gauss("legendre", a , b)
306
307    a = [ 0.06943184,  0.33000948,  0.66999052,  0.93056816]
308    b = [ 0.17392742,  0.32607258,  0.32607258,  0.17392742]
309    run_gauss("legendre01", a , b)
310
311    a = [-0.70710678,  0.70710678]
312    b = [ 0.88622693,  0.88622693]
313    run_gauss("hermite", a , b)
314
315    a = [ -2.02018287,  -0.958572465,   0,           0.958572465,   2.02018287]
316    b = [  0.01995324,   0.39361932,    0.94530872,  0.39361932,    0.01995324]
317    run_gauss("hermite", a , b)
318
319    a = [ 0.41577456,  2.29428036,  6.28994508]
320    b = [ 0.71109301,  0.27851773,  0.01038926]
321    run_gauss("laguerre", a , b)
322
323def test_gauss_quadrature_dynamic(verbose = False):
324    n = 5
325
326    A = mp.randmatrix(2 * n, 1)
327
328    def F(x):
329        r = 0
330        for i in xrange(len(A) - 1, -1, -1):
331            r = r * x + A[i]
332        return r
333
334    def run(qtype, FW, R, alpha = 0, beta = 0):
335        X, W = mp.gauss_quadrature(n, qtype, alpha = alpha, beta = beta)
336
337        a = 0
338        for i in xrange(len(X)):
339            a += W[i] * F(X[i])
340
341        b = mp.quad(lambda x: FW(x) * F(x), R)
342
343        c = mp.fabs(a - b)
344
345        if verbose:
346            print(qtype, c, a, b)
347
348        assert c < 1e-5
349
350    run("legendre", lambda x: 1, [-1, 1])
351    run("legendre01", lambda x: 1, [0, 1])
352    run("hermite", lambda x: mp.exp(-x*x), [-mp.inf, mp.inf])
353    run("laguerre", lambda x: mp.exp(-x), [0, mp.inf])
354    run("glaguerre", lambda x: mp.sqrt(x)*mp.exp(-x), [0, mp.inf], alpha = 1 / mp.mpf(2))
355    run("chebyshev1", lambda x: 1/mp.sqrt(1-x*x), [-1, 1])
356    run("chebyshev2", lambda x: mp.sqrt(1-x*x), [-1, 1])
357    run("jacobi", lambda x: (1-x)**(1/mp.mpf(3)) * (1+x)**(1/mp.mpf(5)), [-1, 1], alpha = 1 / mp.mpf(3), beta = 1 / mp.mpf(5) )
358