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

File trac_9720-random-echelonizable-matrices-v4.patch, 15.3 KB (added by rbeezer, 11 years ago)

Same as v3, but removed routines from global namespace

  • sage/matrix/constructor.py

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