Ticket #9720: trac_9720-random-echelonizable-matrices.patch

File trac_9720-random-echelonizable-matrices.patch, 15.8 KB (added by bwonderly, 11 years ago)
  • sage/matrix/all.py

    # HG changeset patch
    # User Billy Wonderly <bwonderly@ups.edu>
    # Date 1281466684 25200
    # Node ID 3b8d9c6f451f57edc82d0705f19e68fae5b06bd4
    # Parent  426be7b253ad9c3137e733aba26165a3de5b832f
    #9720: add random echelonizable matrices
    
    diff -r 426be7b253ad -r 3b8d9c6f451f sage/matrix/all.py
    a b  
    11from matrix_space import MatrixSpace, is_MatrixSpace
    2 from constructor import matrix, Matrix, random_matrix, diagonal_matrix, identity_matrix, block_matrix, block_diagonal_matrix, jordan_block, zero_matrix
     2from constructor import matrix, Matrix, random_matrix, diagonal_matrix, identity_matrix, block_matrix, block_diagonal_matrix, jordan_block, zero_matrix, random_rref_matrix, random_echelonizable_matrix
    33from matrix import is_Matrix
    44from berlekamp_massey import berlekamp_massey
    55
  • sage/matrix/constructor.py

    diff -r 426be7b253ad -r 3b8d9c6f451f sage/matrix/constructor.py
    a b  
    2525from sage.rings.real_double import RDF
    2626from sage.rings.complex_double import CDF
    2727from sage.rings.integer_ring import ZZ
     28from sage.rings.rational_field import QQ
    2829from sage.misc.misc_c import running_total
    2930from matrix import is_Matrix
    3031from copy import copy
     
    12171218        block[i,i+1]=1
    12181219    return block
    12191220
     1221def random_rref_matrix(num_row, num_col, num_pivots,ring=QQ):
     1222    """
     1223    Generate a matrix already in reduced row-echelon form with the desired row dimension, column dimension, and rank over the designated ring.
     1224
     1225    INPUT:
     1226
     1227    - ``num_row`` - The number of rows desired for the return matrix.
     1228
     1229    - ``num_col`` - The number of columns desired for the return matrix.
     1230
     1231    - ``num_pivots`` - The desired rank of the return matrix.
     1232
     1233    - ``ring`` - The desired ring for the return matrix to be built over.
     1234
     1235    OUTPUT:
     1236
     1237    A matrix in reduced row echelon form with dimensions ``num_row`` x ``num_col`` and a rank of
     1238    ``num_pivot``.
     1239
     1240    EXAMPLES:
     1241
     1242    Matrices generated are in reduced row-echelon form with the default ring of QQ and specified rank. ::
     1243
     1244        sage: A=random_rref_matrix(5,6,4); A # random
     1245        [ 1  0  0 -6  0 -3]
     1246        [ 0  1  0  2  0  3]
     1247        [ 0  0  1 -4  0 -2]
     1248        [ 0  0  0  0  1  3]
     1249        [ 0  0  0  0  0  0]
     1250        sage: A.rank()
     1251        4
     1252        sage: A.base_ring()
     1253        Rational Field
     1254        sage: A==A.rref()
     1255        True
     1256
     1257    Matrices can be generated over other fields. ::
     1258
     1259        sage: B=random_rref_matrix(4,4,3,GF(7)); B # random
     1260        [1 0 0 0]
     1261        [0 1 0 6]
     1262        [0 0 1 4]
     1263        [0 0 0 0]
     1264        sage: B.base_ring()
     1265        Finite Field of size 7
     1266        sage: B==B.rref()
     1267        True
     1268
     1269    TESTS:
     1270
     1271    Row dimension input for a matrix must be whole valued. ::
     1272
     1273       sage: random_rref_matrix(66.9,19,12)
     1274       Traceback (most recent call last):
     1275       ...
     1276       TypeError: inputs for matrix dimensions and rank must be integers.
     1277
     1278    Column dimension input for a matrix must be whole valued. ::
     1279
     1280       sage: random_rref_matrix(31,74.2,6)
     1281       Traceback (most recent call last):
     1282       ...
     1283       TypeError: inputs for matrix dimensions and rank must be integers.
     1284
     1285    Rank input of a matrix must be whole valued. ::
     1286
     1287        sage: random_rref_matrix(120,56,61/2)
     1288        Traceback (most recent call last):
     1289        ...
     1290        TypeError: inputs for matrix dimensions and rank must be integers.
     1291
     1292    Matrices must be generated over exact fields. ::
     1293
     1294       sage: random_rref_matrix(40,88,39,RR)
     1295       Traceback (most recent call last):
     1296       ...
     1297       TypeError: input ring must be exact.
     1298
     1299    Matrices must have the number of pivot columns be less than or equal to the number of rows. ::
     1300
     1301        sage: C=random_rref_matrix(6,4,7); C
     1302        Traceback (most recent call last):
     1303        ...
     1304        ValueError: number of pivot column cannot exceed the number of rows or columns.
     1305
     1306    Matrices must have the number of pivot columns be less than or equal to the number of columns. ::
     1307
     1308        sage: D=random_rref_matrix(1,3,5); D
     1309        Traceback (most recent call last):
     1310        ...
     1311        ValueError: number of pivot column cannot exceed the number of rows or columns.
     1312
     1313    Matrices must have the number of pivot columns be greater than zero. ::
     1314
     1315        sage: random_rref_matrix(5,4,0)
     1316        Traceback (most recent call last):
     1317        ...
     1318        ValueError: matrices must have dimensions and rank greater than zero.
     1319
     1320    AUTHOR:
     1321
     1322    Billy Wonderly (2010-07)
     1323    """
     1324
     1325    import sage.gsl.probability_distribution as pd
     1326    from sage.misc.prandom import randint
     1327
     1328    try:
     1329        num_row=ZZ(num_row)
     1330        num_col=ZZ(num_col)
     1331        num_pivots=ZZ(num_pivots)
     1332    except TypeError:
     1333        raise TypeError("inputs for matrix dimensions and rank must be integers.")
     1334    if not ring.is_exact():
     1335        raise TypeError("input ring must be exact.")
     1336    if num_pivots>num_row or num_pivots>num_col:
     1337        raise ValueError("number of pivot column cannot exceed the number of rows or columns.")
     1338    if num_pivots<=0:
     1339        raise ValueError("matrices must have dimensions and rank greater than zero.")
     1340    else:
     1341        one=ring.one()
     1342        # Create a matrix of the desired size to be modified and then returned.
     1343        return_matrix=matrix(ring,num_row,num_col)
     1344        pivots=[0] #Force first column to be a pivot.
     1345        # Probability distribution for the placement of leading one's.
     1346        pivot_generator=pd.RealDistribution("beta",[1.6,4.3])
     1347        while len(pivots)<num_pivots:
     1348            pivot_column=int(pivot_generator.get_random_element()*num_col)
     1349            if pivot_column not in pivots:
     1350                pivots.append(pivot_column)
     1351        pivots.sort()
     1352        pivot_row=0
     1353        # Use the list of pivot columns to set the pivot entries of the return_matrix to leading ones.
     1354        while pivot_row<num_pivots:
     1355            return_matrix[pivot_row,pivots[pivot_row]]=one
     1356            pivot_row+=1
     1357        if ring==QQ or ring==ZZ:
     1358            # Keep track of the non-pivot columns by using the pivot_index, start at the first column to
     1359            # the right of the initial pivot column, go until the first column to the left of the next
     1360            # pivot column.
     1361            for pivot_index in range(num_pivots-1):
     1362                for non_pivot_column_index in range(pivots[pivot_index]+1,pivots[pivot_index+1]):
     1363                    entry_generator1=pd.RealDistribution("beta",[6,4])
     1364                    # Experimental distribution used to generate the values.
     1365                    for non_pivot_column_entry in range(pivot_index+1):
     1366                        sign1=(2*randint(0,1)-1)
     1367                        return_matrix[non_pivot_column_entry,non_pivot_column_index]=sign1*int(entry_generator1.get_random_element()*((1-non_pivot_column_entry/return_matrix.ncols())*7))
     1368            # Use index to fill entries of the columns to the right of the last pivot column.
     1369            for rest_non_pivot_column in range(pivots[num_pivots-1]+1,num_col):
     1370                entry_generator2=pd.RealDistribution("beta",[2.6,4])
     1371                # experimental distribution to generate small values.
     1372                for rest_entries in range(num_pivots):
     1373                    sign2=(2*randint(0,1)-1)
     1374                    return_matrix[rest_entries,rest_non_pivot_column]=sign2*int(entry_generator2.get_random_element()*5)
     1375        else:
     1376            for pivot_index in range(num_pivots-1):
     1377                for non_pivot_column_index in range(pivots[pivot_index]+1,pivots[pivot_index+1]):
     1378                    for non_pivot_column_entry in range(pivot_index+1):
     1379                            return_matrix[non_pivot_column_entry,non_pivot_column_index]=ring.random_element()
     1380            for rest_non_pivot_column in range(pivots[num_pivots-1]+1,num_col):
     1381                for rest_entries in range(num_pivots):
     1382                    return_matrix[rest_entries,rest_non_pivot_column]=ring.random_element()
     1383    return return_matrix
     1384
     1385def random_echelonizable_matrix(rows,columns,rank,upper_bound=None,ring=QQ):
     1386    """
     1387    Generate a matrix of a desired size and rank, over a desired field, whose reduced
     1388    row-echelon form has only whole values.
     1389
     1390    INPUT:
     1391
     1392    - ``rows`` - The row dimension of the desired matrix, taking only integers greater than zero.
     1393
     1394    - ``columns`` - The column dimension of the desired matrix, taking only integers greater than zero.
     1395
     1396    - ``rank`` - The desired rank, or number of leading ones, for the return matrix.
     1397
     1398    - ``upper_bound`` - If designated, size control of the matrix entries is desired, ``upper_bound`` is 1 + the maximum value entries can achieve. 
     1399      If None, no size control occurs. (default: None)
     1400
     1401    - ``ring`` - The desired ring for the matrix to be generated over (default: QQ).
     1402
     1403
     1404    OUTPUT:
     1405
     1406    A matrix not in reduced row-echelon form with the desired dimensions and properties.
     1407
     1408    EXAMPLES:
     1409
     1410    Generated matrices have the desired dimensions, rank and entry size. The matrix in reduced row-echelon form has only integer entries. ::
     1411
     1412        sage: A=random_echelonizable_matrix(5,6,4,upper_bound=40,ring=QQ); A # random
     1413        [  1  -1   1  -3  -4   6]
     1414        [  5  -4   0   8   4  19]
     1415        [ -3   3  -2   4   7 -16]
     1416        [ -4   5  -7  26  31 -31]
     1417        [  2  -3   4 -11 -14  17]
     1418        sage: A.rank()
     1419        4
     1420        sage: max(map(abs,A.list()))<40
     1421        True
     1422        sage: A.rref()==A.rref().change_ring(ZZ)
     1423        True
     1424
     1425    An example with default settings (ring=QQ, no entry size control). ::
     1426
     1427        sage: C=random_echelonizable_matrix(6,7,5); C # random
     1428        [  1   0   5  -2 -26 -16   0]
     1429        [ -3   1 -19   6  97  61   1]
     1430        [  0   4 -15  -1  71  50   3]
     1431        [  2   4  -9   0  39  25   8]
     1432        [  2   2   3  -3 -18  -9   3]
     1433        [ -3  -4  -2  14  14  -6   4]
     1434        sage: C.rank()
     1435        5
     1436        sage: C.rref()==C.rref().change_ring(ZZ)
     1437        True
     1438
     1439    A matrix without size control may have very large entry sizes. ::
     1440
     1441        sage: D=random_echelonizable_matrix(7,8,6,ring=ZZ); D # random
     1442        [    9   -53  -255    45 -1519  4043  9819  3324]
     1443        [    3   -14   -64     8  -369   972  2350   810]
     1444        [    2   -14   -65     9  -377  1000  2420   829]
     1445        [    4   -24  -116    21  -693  1846  4485  1516]
     1446        [   -3    14    68   -16   426 -1134 -2767  -919]
     1447        [   -5    21    92   -13   548 -1432 -3466 -1183]
     1448        [    1    -9   -42     7  -254   670  1624   547]
     1449
     1450    Matrices can be generated over a different ring. ::
     1451
     1452        sage: F.<a>=GF(2^3)
     1453        sage: B=random_echelonizable_matrix(4,5,4,None,F); B # random
     1454        [      a + 1 a^2 + a + 1         a^2           0           1]
     1455        [          1           a       a + 1         a^2     a^2 + a]
     1456        [          a         a^2 a^2 + a + 1       a + 1           1]
     1457        [          1           0 a^2 + a + 1           0     a^2 + 1]
     1458        sage: B.rank()
     1459        4
     1460
     1461    Square matrices over ZZ or QQ with full rank are unimodular. ::
     1462
     1463        sage: E=random_echelonizable_matrix(7,7,7); E # random
     1464        [  1  -1   5  12 -24 -41  47]
     1465        [  0   1  -1   3   0 -11  40]
     1466        [  1  -1   6   6 -19 -20 -11]
     1467        [ -2   1 -10 -12  35  44  -4]
     1468        [  3  -1   9   7 -35 -40 -18]
     1469        [  0   0   0  -4   4  13 -32]
     1470        [  3  -3  11   6 -33 -31 -35]
     1471        sage: det(E)
     1472        1
     1473
     1474    TESTS:
     1475
     1476    Matrices must have a row dimension greater than zero. ::
     1477
     1478        sage: random_echelonizable_matrix(0,5,3)
     1479        Traceback (most recent call last):
     1480        ...
     1481        ValueError: matrices must have dimensions and rank greater than zero.
     1482
     1483    Matrices must have a column dimension greater than zero. ::
     1484   
     1485        sage: random_echelonizable_matrix(9,0,4)
     1486        Traceback (most recent call last):
     1487        ...
     1488        ValueError: matrices must have dimensions and rank greater than zero.
     1489
     1490    Matrices must have a rank greater than zero. ::
     1491
     1492        sage: random_echelonizable_matrix(3,4,0)
     1493        Traceback (most recent call last):
     1494        ...
     1495        ValueError: matrices must have dimensions and rank greater than zero.
     1496
     1497    AUTHOR:
     1498
     1499    Billy Wonderly (2010-07)
     1500    """
     1501
     1502    from sage.misc.prandom import randint
     1503
     1504    if rows<=0 or columns<=0 or rank<=0:
     1505        raise ValueError("matrices must have dimensions and rank greater than zero.")
     1506    # Entries of matrices over the ZZ or QQ can get large, entry size is regulated by finding the largest
     1507    # entry of the resultant matrix after addition of scalar multiple of a row.
     1508    if ring==QQ or ring==ZZ:
     1509        matrix=random_rref_matrix(rows,columns,rank,ring)
     1510        matrix_copy=copy(matrix)
     1511        # If upper_bound is not set, don't control entry size.
     1512        if upper_bound==None:
     1513            upper_bound=50
     1514            size=False
     1515        else:
     1516            size=True
     1517        if size:
     1518            for pivots in range(len(matrix.pivots())-1,-1,-1):
     1519            # keep track of the pivot column positions from the pivot column with the largest index to
     1520            # the one with the smallest.
     1521                row_index=0
     1522                while row_index<rows:
     1523                    # To each row in a pivot column add a scalar multiple of the pivot row.
     1524                    if pivots!=row_index:
     1525                    # To ensure a leading one is not removed by the addition of the pivot row by its
     1526                    # additive inverse.
     1527                        matrix_copy=matrix.with_added_multiple_of_row(row_index,matrix.pivot_rows()[pivots],randint(-5,5))
     1528                        # Range for scalar multiples determined experimentally.
     1529                    if max(map(abs,matrix_copy.list()))<upper_bound:
     1530                    # Continue if the the largest entry after a row operation is within the bound.
     1531                        matrix=matrix_copy
     1532                        row_index+=1
     1533            # The leading one in row one has not been altered, so add a scalar multiple of a random row
     1534            # to row one.
     1535            row1=0
     1536            if rows>1:
     1537                while row1<1:
     1538                    matrix_copy=matrix.with_added_multiple_of_row(0,randint(1,rows-1),randint(-3,3))
     1539                    if max(map(abs,matrix_copy.list()))<upper_bound:
     1540                        matrix=matrix_copy
     1541                        row1+=1
     1542        # If size control is not desired, the routine will run slightly faster, particularly with large matrices.
     1543        else:
     1544            for pivots in range(rank-1,-1,-1):
     1545                row_index=0
     1546                while row_index<rows:
     1547                    if pivots==row_index:
     1548                        row_index+=1
     1549                    if pivots!=row_index and row_index!=rows:
     1550                        matrix.add_multiple_of_row(row_index,matrix.pivot_rows()[pivots],randint(-5,5))
     1551                        row_index+=1
     1552            if rows>1:
     1553                matrix.add_multiple_of_row(0,randint(1,rows-1),randint(-3,3))
     1554    # If the matrix generated over a different ring, random elements from the designated ring are used as and
     1555    # the routine is run similarly to the size unchecked version for rationals and integers.
     1556    else:
     1557        matrix=random_rref_matrix(rows,columns,rank,ring)
     1558        matrix_copy=copy(matrix)
     1559        for pivots in range(rank-1,-1,-1):
     1560            row_index=0
     1561            while row_index<rows:
     1562                if pivots==row_index:
     1563                    row_index+=1
     1564                if pivots!=row_index and row_index!=rows:
     1565                    matrix.add_multiple_of_row(row_index,matrix.pivot_rows()[pivots],ring.random_element())
     1566                    row_index+=1
     1567        if rows>1:
     1568            matrix.add_multiple_of_row(0,randint(1,rows-1),ring.random_element())
     1569    return matrix