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

File trac_9720-random-echelonizable-matrices-v3.patch, 15.8 KB (added by bwonderly, 11 years ago)

stand-alone version 3 of the same patch, with import statement change and documentation revision

  • sage/matrix/all.py

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