This file is indexed.

/usr/share/pyshared/cogent/struct/knots.py is in python-cogent 1.5.3-2.

This file is owned by root:root, with mode 0o644.

The actual contents of the file can be viewed below.

   1
   2
   3
   4
   5
   6
   7
   8
   9
  10
  11
  12
  13
  14
  15
  16
  17
  18
  19
  20
  21
  22
  23
  24
  25
  26
  27
  28
  29
  30
  31
  32
  33
  34
  35
  36
  37
  38
  39
  40
  41
  42
  43
  44
  45
  46
  47
  48
  49
  50
  51
  52
  53
  54
  55
  56
  57
  58
  59
  60
  61
  62
  63
  64
  65
  66
  67
  68
  69
  70
  71
  72
  73
  74
  75
  76
  77
  78
  79
  80
  81
  82
  83
  84
  85
  86
  87
  88
  89
  90
  91
  92
  93
  94
  95
  96
  97
  98
  99
 100
 101
 102
 103
 104
 105
 106
 107
 108
 109
 110
 111
 112
 113
 114
 115
 116
 117
 118
 119
 120
 121
 122
 123
 124
 125
 126
 127
 128
 129
 130
 131
 132
 133
 134
 135
 136
 137
 138
 139
 140
 141
 142
 143
 144
 145
 146
 147
 148
 149
 150
 151
 152
 153
 154
 155
 156
 157
 158
 159
 160
 161
 162
 163
 164
 165
 166
 167
 168
 169
 170
 171
 172
 173
 174
 175
 176
 177
 178
 179
 180
 181
 182
 183
 184
 185
 186
 187
 188
 189
 190
 191
 192
 193
 194
 195
 196
 197
 198
 199
 200
 201
 202
 203
 204
 205
 206
 207
 208
 209
 210
 211
 212
 213
 214
 215
 216
 217
 218
 219
 220
 221
 222
 223
 224
 225
 226
 227
 228
 229
 230
 231
 232
 233
 234
 235
 236
 237
 238
 239
 240
 241
 242
 243
 244
 245
 246
 247
 248
 249
 250
 251
 252
 253
 254
 255
 256
 257
 258
 259
 260
 261
 262
 263
 264
 265
 266
 267
 268
 269
 270
 271
 272
 273
 274
 275
 276
 277
 278
 279
 280
 281
 282
 283
 284
 285
 286
 287
 288
 289
 290
 291
 292
 293
 294
 295
 296
 297
 298
 299
 300
 301
 302
 303
 304
 305
 306
 307
 308
 309
 310
 311
 312
 313
 314
 315
 316
 317
 318
 319
 320
 321
 322
 323
 324
 325
 326
 327
 328
 329
 330
 331
 332
 333
 334
 335
 336
 337
 338
 339
 340
 341
 342
 343
 344
 345
 346
 347
 348
 349
 350
 351
 352
 353
 354
 355
 356
 357
 358
 359
 360
 361
 362
 363
 364
 365
 366
 367
 368
 369
 370
 371
 372
 373
 374
 375
 376
 377
 378
 379
 380
 381
 382
 383
 384
 385
 386
 387
 388
 389
 390
 391
 392
 393
 394
 395
 396
 397
 398
 399
 400
 401
 402
 403
 404
 405
 406
 407
 408
 409
 410
 411
 412
 413
 414
 415
 416
 417
 418
 419
 420
 421
 422
 423
 424
 425
 426
 427
 428
 429
 430
 431
 432
 433
 434
 435
 436
 437
 438
 439
 440
 441
 442
 443
 444
 445
 446
 447
 448
 449
 450
 451
 452
 453
 454
 455
 456
 457
 458
 459
 460
 461
 462
 463
 464
 465
 466
 467
 468
 469
 470
 471
 472
 473
 474
 475
 476
 477
 478
 479
 480
 481
 482
 483
 484
 485
 486
 487
 488
 489
 490
 491
 492
 493
 494
 495
 496
 497
 498
 499
 500
 501
 502
 503
 504
 505
 506
 507
 508
 509
 510
 511
 512
 513
 514
 515
 516
 517
 518
 519
 520
 521
 522
 523
 524
 525
 526
 527
 528
 529
 530
 531
 532
 533
 534
 535
 536
 537
 538
 539
 540
 541
 542
 543
 544
 545
 546
 547
 548
 549
 550
 551
 552
 553
 554
 555
 556
 557
 558
 559
 560
 561
 562
 563
 564
 565
 566
 567
 568
 569
 570
 571
 572
 573
 574
 575
 576
 577
 578
 579
 580
 581
 582
 583
 584
 585
 586
 587
 588
 589
 590
 591
 592
 593
 594
 595
 596
 597
 598
 599
 600
 601
 602
 603
 604
 605
 606
 607
 608
 609
 610
 611
 612
 613
 614
 615
 616
 617
 618
 619
 620
 621
 622
 623
 624
 625
 626
 627
 628
 629
 630
 631
 632
 633
 634
 635
 636
 637
 638
 639
 640
 641
 642
 643
 644
 645
 646
 647
 648
 649
 650
 651
 652
 653
 654
 655
 656
 657
 658
 659
 660
 661
 662
 663
 664
 665
 666
 667
 668
 669
 670
 671
 672
 673
 674
 675
 676
 677
 678
 679
 680
 681
 682
 683
 684
 685
 686
 687
 688
 689
 690
 691
 692
 693
 694
 695
 696
 697
 698
 699
 700
 701
 702
 703
 704
 705
 706
 707
 708
 709
 710
 711
 712
 713
 714
 715
 716
 717
 718
 719
 720
 721
 722
 723
 724
 725
 726
 727
 728
 729
 730
 731
 732
 733
 734
 735
 736
 737
 738
 739
 740
 741
 742
 743
 744
 745
 746
 747
 748
 749
 750
 751
 752
 753
 754
 755
 756
 757
 758
 759
 760
 761
 762
 763
 764
 765
 766
 767
 768
 769
 770
 771
 772
 773
 774
 775
 776
 777
 778
 779
 780
 781
 782
 783
 784
 785
 786
 787
 788
 789
 790
 791
 792
 793
 794
 795
 796
 797
 798
 799
 800
 801
 802
 803
 804
 805
 806
 807
 808
 809
 810
 811
 812
 813
 814
 815
 816
 817
 818
 819
 820
 821
 822
 823
 824
 825
 826
 827
 828
 829
 830
 831
 832
 833
 834
 835
 836
 837
 838
 839
 840
 841
 842
 843
 844
 845
 846
 847
 848
 849
 850
 851
 852
 853
 854
 855
 856
 857
 858
 859
 860
 861
 862
 863
 864
 865
 866
 867
 868
 869
 870
 871
 872
 873
 874
 875
 876
 877
 878
 879
 880
 881
 882
 883
 884
 885
 886
 887
 888
 889
 890
 891
 892
 893
 894
 895
 896
 897
 898
 899
 900
 901
 902
 903
 904
 905
 906
 907
 908
 909
 910
 911
 912
 913
 914
 915
 916
 917
 918
 919
 920
 921
 922
 923
 924
 925
 926
 927
 928
 929
 930
 931
 932
 933
 934
 935
 936
 937
 938
 939
 940
 941
 942
 943
 944
 945
 946
 947
 948
 949
 950
 951
 952
 953
 954
 955
 956
 957
 958
 959
 960
 961
 962
 963
 964
 965
 966
 967
 968
 969
 970
 971
 972
 973
 974
 975
 976
 977
 978
 979
 980
 981
 982
 983
 984
 985
 986
 987
 988
 989
 990
 991
 992
 993
 994
 995
 996
 997
 998
 999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217
1218
1219
1220
1221
1222
1223
1224
1225
1226
1227
1228
1229
1230
1231
1232
1233
1234
1235
1236
1237
1238
1239
1240
1241
1242
1243
1244
1245
1246
1247
1248
1249
1250
1251
1252
1253
1254
1255
1256
1257
1258
1259
1260
1261
1262
1263
1264
1265
1266
1267
1268
1269
1270
1271
1272
1273
1274
1275
1276
1277
1278
1279
1280
1281
1282
1283
1284
1285
1286
1287
1288
1289
1290
1291
1292
1293
1294
1295
1296
1297
1298
1299
1300
1301
1302
1303
1304
1305
1306
1307
1308
1309
1310
1311
1312
1313
1314
1315
1316
1317
1318
1319
1320
1321
1322
1323
1324
1325
1326
1327
1328
1329
1330
1331
1332
1333
1334
1335
1336
1337
1338
1339
1340
1341
1342
1343
1344
1345
1346
1347
1348
1349
1350
1351
1352
1353
1354
1355
1356
1357
1358
1359
1360
1361
1362
1363
1364
1365
1366
1367
1368
1369
1370
1371
1372
1373
1374
1375
1376
1377
1378
1379
1380
1381
1382
1383
1384
1385
1386
1387
1388
1389
1390
1391
1392
1393
1394
1395
1396
1397
1398
1399
1400
1401
1402
1403
1404
1405
1406
1407
1408
1409
1410
1411
1412
1413
1414
1415
1416
1417
1418
1419
1420
1421
1422
1423
1424
1425
1426
1427
1428
1429
1430
1431
1432
1433
1434
1435
1436
1437
1438
1439
1440
1441
1442
1443
1444
1445
1446
1447
1448
1449
1450
1451
1452
1453
1454
1455
1456
1457
1458
1459
1460
1461
1462
1463
1464
1465
1466
1467
1468
1469
1470
1471
1472
1473
1474
1475
1476
1477
1478
1479
1480
1481
1482
1483
1484
1485
1486
1487
1488
1489
1490
1491
1492
1493
1494
1495
1496
1497
1498
1499
1500
1501
1502
1503
1504
1505
1506
1507
1508
1509
1510
1511
1512
1513
1514
1515
1516
1517
1518
1519
1520
1521
1522
1523
1524
1525
1526
1527
1528
1529
1530
1531
1532
1533
1534
1535
1536
1537
1538
1539
1540
1541
1542
1543
1544
1545
1546
1547
1548
1549
1550
1551
1552
1553
1554
1555
1556
1557
1558
1559
1560
1561
1562
1563
1564
1565
1566
1567
1568
1569
1570
1571
1572
1573
1574
1575
1576
1577
1578
1579
1580
1581
1582
1583
1584
1585
1586
1587
1588
1589
1590
1591
1592
1593
1594
1595
1596
1597
1598
1599
1600
1601
1602
1603
1604
1605
1606
1607
1608
1609
1610
1611
1612
1613
1614
1615
1616
1617
1618
1619
1620
1621
1622
1623
1624
1625
1626
1627
1628
1629
1630
1631
1632
1633
1634
1635
1636
1637
1638
1639
1640
1641
1642
1643
1644
1645
1646
1647
1648
1649
1650
1651
1652
1653
1654
1655
1656
1657
1658
1659
1660
1661
1662
1663
1664
1665
1666
1667
1668
1669
1670
1671
1672
1673
1674
1675
1676
1677
1678
1679
1680
1681
1682
1683
1684
1685
1686
1687
1688
1689
1690
1691
1692
1693
1694
1695
1696
1697
1698
1699
1700
1701
#!/usr/bin/env python
# knots.py

"""Contains code related to RNA (secondary) structure and pseudoknots.

Specifically, this module contains several methods to remove pseudoknots
from RNA structures. Pseudoknot removal is discussed in the following paper:

S. Smit, K. Rother, J. Heringa, and R. Knight
Manuscript in preparation.

If you use this code in your work, please cite this publication (in addition
to the PyCogent publication). If you need to cite the paper before submission,
please contact the author of this module.

Six functions are provided (see documentation of each function for
a more detailed description):

opt_all -- optimization approach that calculates all nested structures
    that optimize some value (e.g. keep the maximum number of base pairs)
conflict_elimination -- Removes pseudoknots from a structure by eliminating
    conflicting base pairs one by one. Two functions to determine which 
    paired region should be removed next are provided: max_conlficts
    and min_gain.
inc_order -- creates a nested structure by adding non-conflicting paired
    regions one by one to the solution; paired regions are processed
    from 5' to 3' start point or from 3' to 5' end point.
inc_length -- creates a nested structure by adding non-conflicting paired
    regions one at the time, starting with the longest region working towards
    the shortest region.
inc_range -- generates a nested structure by adding non-conflicting paired
    retions one at the time, starting with short-range interactions working
    towards long-range interactions.

These six functions represent the core objective of this module.
Two convenience functions supporting the opt_all function are added:
opt_single_random, and opt_single_property
There is also a modified version of the original Nussinov-Jacobson
algorithm present which is restricted to the given list of base pairs:
nussinov_restricted

In addition, the following supporting objects and functions are present:

PairedRegion -- object that represents a paired region in an RNA structure.
    A paired region is an uninterrupted stretch of base pairs with positions
    [(i,j),(i+1, j-1),(i+2,j-2), ...].
PairedRegions -- object (basically a list) that stores a collection of 
    PairedRegion objects. This is an alternative way of representing an
    RNA structure, where basically stretches of base pairs are condensed
    into PairedRegion objects.
PairedRegionFromPairs -- Factory function to create a PairedRegion object
    from a Pairs object (cogent.struct.rna2d)
PairedRegionsFromPairs -- Factory function to create a PairedRegions object
    from a Pairs object (cogent.struct.rna2d)
ConflictMatrix -- object to store a matrix of conflicts between different 
    paired regions. Row and Column indices correspond to PairedRegion IDs, 
    values in the matrix are True (if the regions conflict), and False (if
    the regions don't conflict).

Other smaller helper functions are:
contains_true, empty_matrix, pick_multi_best, dp_matrix_multi,
matrix_solutions, and add_back_non_conflicting
See their docstrings for detailed documentation.

NOTE:
None of the methods provided can handle overlapping base pairs (i.e. base
A pairs with B, and A pairs with C), because in that case it is unclear
whether the first or second pair should be kept (also, you could get 
different results based on the order of the base pairs). In general,
one removes pseudoknots in order to represent the structure in dot-bracket
format. If the list contains conflict, this is impossible anyway, so removing
the pseudoknots would not help. It is the responsibility of the user to 
remove overlapping base pairs before trying to obtain a nested structure.

"""

from __future__ import division
from random import choice
from numpy import sum, average, zeros
from cogent.struct.rna2d import Pairs
from cogent.util.dict2d import Dict2D

__author__ = "Sandra Smit"
__copyright__ = "Copyright 2007-2012, The Cogent Project"
__contributors__ = ["Sandra Smit, Rob Knight"]
__license__ = "GPL"
__version__ = "1.5.3"
__maintainer__ = "Sandra Smit"
__email__ = "sandra.smit@colorado.edu"
__status__ = "Production"


class PairedRegion(object):
    """Store an uninterrupted list of base pairs by start, end, and length
    
    A paired region (a.k.a. ladder or (helical) region) is a stretch of
        perfectly nested base pairs with positions:
        [(m,n), (m+1,n-1), (m+2,n-2),...]
    This object is very similar to the Stem object in cogent.struct.rna2d.
        In addition to the start, end, and length it stores the actual 
        base pairs, and it stores a region ID. It has many more methods
        than the Stem object.
    This object performs no error checking. You can specify and End 
        before a Start point, or a Start and End which are closer 
        to each other than 2 times the Length.
    """

    def __init__(self, Start, End, Length, Id=None):
        """Initialize a new PairedRegion object

        Start -- int, specifying the starting index of the paired region
            in the sequence. This is the 5' side of the 5' halfregion.
        End -- int, specifying the end index of the paired region
            in the sequence. This is the 3' side of the 3' halfregion.
        Length -- int, specifying the length of the paired region, i.e. the
            number of base pairs in the region. 
        Id -- string or int, unique identifier for this PairedRegion.

        During initialization a Pairs object is created. The first pair is
        always (Start,End), additional pairs add up from the Start point and
        down from the End point, so (Start+1,End-1), (Start+2, End-2) etc.
        """
        self.Start = Start
        self.End = End
        if Length < 1:
            raise ValueError(\
                "PairedRegion should contain at least one base pair")
        self.Length = Length
        self.Id = Id
        self.Pairs = Pairs()
        for i in range(self.Length):
            self.Pairs.append((self.Start+i, self.End-i)) 

    def __str__(self):
        """Return string representation of PairedRegion, list of Pairs
        """
        return str(self.Pairs)

    def __len__(self):
        """Return Length of the PairedRegion, i.e. the number of base pairs
        """
        return self.Length

    def __eq__(self, other):
        """Compares Pairs and IDs
       
        other -- PairedRegion object

        If IDs are not set (both None), Pairs is the only criterion. 
        """
        if self.Pairs == other.Pairs and self.Id == other.Id:
            return True
        return False

    def __ne__(self, other):
        """Return True if two PairedRegion objects differ

        other -- PairedRegion object
        """
        return not self == other

    def upstream(self):
        """Return list of upstream positions in self from 5' to 3'
        """
        return [i for (i,j) in self.Pairs]

    def downstream(self):
        """Return list of downstream positions in self from 5' to 3'
        """
        result = [j for (i,j) in self.Pairs]
        result.reverse()
        return result

    def paired(self):
        """Return sorted list of paired positions in this region
        """
        result = self.upstream() + self.downstream()
        result.sort()
        return result

    def range(self):
        """Return the range of this region
        
        The range of the region is the number of bases between the
        highest upstream and the lowest downstream position.
        (i.e. the number of unpaired bases in the hairpin if this
        were the only paired region in the structure)
        For example: the range of region with start=3, end=10, and len=2
            would be 4.
        Performs no error checking. If region overlaps with itself,
            a negative number will be returned...
        """
        return min(self.downstream()) - max(self.upstream()) - 1

    def overlapping(self, other):
        """Returns True if two regions overlap

        other -- PairedRegion object

        Two regions overlap if there is at least one base which
            is a member of both regions. (Definition from Studnicka 1978)
        Identical regions are overlapping.
        """
        ref_pos = dict.fromkeys(self.paired())
        for pos in other.paired():
            if pos in ref_pos:
                return True
        return False

    def conflicting(self, other):
        """Return True if the regions are conflicting, False if they are nested

        other -- PairedRegion object

        Two paired regions are conflicting if they are organized in a 
            knotted fashion. This means both other.Start and other.End 
            have to be between self.Start and self.End or both shouldn't be.
            See for example Studnicka 1978 for defintion, or any other
            paper with a general pseudoknot definiton in there.
        Overlapping regions cause an error, because you can't determine
            whether they are conflicting or not. However, two identical
            regions are defined as NOT conflicting, even though they are
            overlapping, and thus False is returned. For non-identical,
            but overlapping regions an error will be raised.
        """
        if self == other: # equal blocks
            return False
        # if not equal, but overlapping, raise error
        if self.overlapping(other):
            raise ValueError("Can only handle non-overlapping regions")
        
        if (other.Start > self.End and other.End > self.End) or\
            (self.Start > other.End and self.End > other.End):
            return False
        if (self.Start < other.Start < self.End and\
            self.Start < other.End < self.End) or\
            (other.Start < self.Start < other.End and\
            other.Start < self.End < other.End):
            return False
        return True

    def score(self, scoring_function):
        """Sets self.Score to value of scoring function applies to self

        scoring_function -- function that can be applied to a PairedRegion
            object and returns a numerical value

        Note: this method has no return value, it sets a property of the 
            instance.
        """
        self.Score = scoring_function(self)

def PairedRegionFromPairs(pairs, Id=None):
    """Return new PairedRegion object from Pairs object

    pairs -- Pairs object or list of tuples with up and downstream positions.
    Id -- string or int, unique identifier of this region.

    This is a factory function to create a PairedRegion object
        from a Pairs object. It assumes the pairs are fully nested
        with positions [(m,n), (m+1,n-1), (m+2,n-2),...].
    The function doesn't validate the input pairs, so if the assumtion
        does not hold, the Pairs in the resulting PairedRegion might
        differ from the input pairs.
    It makes the pairs directed and sorts them, then it extracts the 
        Start, End and Length and initialized a new PairedRegion with 
        those parameters.
    """
    if not pairs:
        raise ValueError("PairedRegion should contain at least one pair")
    # preprocess
    raw_pairs = Pairs(pairs)
    if raw_pairs.hasConflicts():
        raise ValueError("Cannot handle pairs with conflicts")
    p = raw_pairs.directed()
    p.sort()
    # initialize variables
    start = p[0][0]
    end = p[0][1]
    length = len(p)
    return PairedRegion(start, end, length, Id=Id) 


class PairedRegions(list):
    """Stores a list of PairedRegion objects
    
    A PairedRegions object is a condensed way of looking at an RNA structure,
        where continuous stretches of base pairs are collapsed into 
        PairedRegion objects. See the documentation on PairedRegion for more
        details.
    """
    
    def __init__(self, regions=None):
        """Initialize new PairedRegions object

        regions -- list of PairedRegion objects

        The object is meant to store a list of PairedRegion objects,
            however it is very light-weight and does not perform any
            validation on the input.
        """
        if regions is None:
            regions = []
        self[:] = list(regions)

    def __str__(self):
        """Return string representation of PairedRegions object

        Each PairedRegion is presented as Id:Start,End,Length;
            A PairedRegions object is presented as '(' + space-delimited
            list of PairedRegion objects + ')'
            For example: (A:2,10,2; B:12,20,3;)
        """
        result = []
        for i in self:
            result.append("%s:%s,%s,%s;"%(i.Id,i.Start,i.End,i.Length))
        return '('+' '.join(result)+')'
    
    def __eq__(self, other):
        """Return True if two PairedRegions objects are equal

        other -- PairedRegions object

        Two regions are equal if they have the same length and contain
            the same PairedRegion objects.
        """
        if len(self) != len(other):
            return False
        for i in self:
            if i not in other:
                return False
        return True

    def __ne__(self, other):
        """Return True if two PairedRegions objects are different

        other -- PairedRegions object
        """
        return not self == other

    def byId(self):
        """Return dict of {ID: PairedRegion}

        This function only works if the IDs for each region in self
            are unique. If multiple regions with the same ID are found, 
            an error is raised.
        """
        result = {}
        for pr in self: # for every paired region
            if pr.Id in result:
                raise ValueError("Duplicate key found")
            result[pr.Id] = pr
        return result

    def numberOfRegions(self):
        """Return the number of PairedRegion objects in the list
        """
        return len(self)

    def totalLength(self):
        """Return the cumulative length of all PairedRegion objects in self
        
        The totalLength is the total number of base pairs in this
            PairedRegions object. So, it adds the number of pairs in each
            PairedRegion in the list.
        """
        if self:
            return sum(map(len, self))
        else:
            return 0

    def totalScore(self):
        """Return sum of Score values of each PairedRegion in self

        This method simply adds all the Score attributes (!= None)
            for each PairedRegion in this PairedRegions object.
        """
        score = 0
        for pr in self:
            try:
                score += pr.Score
            except AttributeError:
                raise ValueError("Score not set for %s"%(str(self)))
            except TypeError:
                raise ValueError("Score should be numerical, but is %s"\
                    %(pr.Score))
        return score

    def toPairs(self):
        """Return Pairs object containing all the pairs in each PairedRegion

        This method does not validate the pairs. It simply adds all the 
            pairs in each PairedRegion to the result. Pairs might occur twice
            in the result. The resulting Pairs object is sorted.
        """
        result = Pairs()
        for pr in self:
            result.extend(pr.Pairs)
        result.sort()
        return result

    def byStartEnd(self):
        """Return dict of {(pr.Start, pr.End): pr}

        Keys in the dictionary are tuples of start and end positions,
            the values are the PairedRegion objects themselves. If a 
            Start/End combination is already in the dictionary, an error is 
            raised.
        """
        result = {}
        for pr in self:
            se = (pr.Start, pr.End)
            if se in result:
                raise ValueError("Duplicate key found: %s"%(str(se)))
            result[se] = pr
        return result
        #return dict([((pr.Start, pr.End), pr) for pr in self])

    def lowestStart(self):
        """Return lowest begin value of any PairedRegion in the list
        
        It lists all the Start values for all the PairedRegion objects
            in the list and returns the lowest value. If there are no
            regions in self, None is returned.
        """
        start_values = [pr.Start for pr in self]
        if not start_values:
            return None
        else:
            return min(start_values)

    def highestEnd(self):
        """Return highest end value of any PairedRegion in the list

        It lists all the End values for all the PairedRegion objects
            in the list and returns the highest value. If there are no
            regions in self, None is returned.
        """
        end_values = [pr.End for pr in self]
        if not end_values:
            return None
        else:
            return max(end_values)

    def sortedIds(self):
        """Return sorted list of region IDs
        """
        all_ids = [pr.Id for pr in self]
        all_ids.sort()
        return all_ids

    def upstream(self):
        """Return sorted list of upstream positions
        """
        result = []
        for pr in self:
            result.extend(pr.upstream())
        result.sort()
        return result

    def downstream(self):
        """Return sorted list of downstream positions
        """
        result = []
        for pr in self:
            result.extend(pr.downstream())
        result.sort()
        return result

    def pairedPos(self):
        """Return sorted list of all paired positions
        """
        result = self.upstream() + self.downstream()
        result.sort()
        return result

    def boundaries(self):
        """Return sorted list of all start and end points
        """
        result = []
        for pr in self:
            result.append(pr.Start)
            result.append(pr.End)
        result.sort()
        return result

    def enumeratedBoundaries(self):
        """Return dict of {boundary_index: boundary}

        Return value is dictionary created from tuples from the enumeration
        of all boundaries.
        """
        return dict(enumerate(self.boundaries()))

    def invertedEnumeratedBoundaries(self):
        """Return dict of {boundary_value: boundary_idx}

        Boundary values are all the start and end points of the paired
            regions in self. Should be unique, otherwise an error is raised.
        Boundary indices are the indices assigned to each start and end point
            during an enumeration of the sorted list. 
        Overall the result is the inverted dictionary of the result of the
            enumeratedBoundaries method.
        """
        eb = self.enumeratedBoundaries()
        result = {}
        for boundary_idx, boundary_value in eb.items():
            if boundary_value in result:
                raise ValueError(\
                    "Boundary value %s is not unique"%(boundary_value))
            result[boundary_value] = boundary_idx
        #result = dict([(v,k) for k,v in eb.items()])
        return result

    def merge(self, other):
        """Merge two PairedRegions objects together

        other -- PairedRegions object

        Duplicate PairedRegion objects are stored only once.
        This methods used the PairedRegion IDs to check for duplications.
        """
        result = PairedRegions()
        seen = {}
        for pr in self+other:
            if pr.Id not in seen:
                result.append(pr)
                seen[pr.Id] = True
        return result

    def conflicting(self, cm=None):
        """Return PairedRegions obj containing regions involved in a conflict

        cm -- ConflictMatrix for this PairedRegions object.

        This method only works if the PairedRegion objects have unique IDs, 
            because a conflict matrix is constructed. This behavior can be 
            changed...
        See PairedRegion.conflicting() for a definition of conflicting
            paired regions.
        """
        if cm is None:
            cm = ConflictMatrix(self)
        id_to_pr = self.byId()

        result = PairedRegions()
        for pr_id in cm.conflicting():
            result.append(id_to_pr[pr_id])
        return result

    def nonConflicting(self, cm=None):
        """Return new PairedRegions object containing non-conflicing regions

        cm -- ConflictMatrix for this PairedRegions object.

        Two PairedRegion objects do not conflict when they are organized
            in a nested fashion.
        """
        if cm is None:
            cm = ConflictMatrix(self)
        id_to_pr = self.byId()

        result = PairedRegions()
        for pr_id in cm.nonConflicting():
            result.append(id_to_pr[pr_id])
        return result

    def conflictCliques(self, cm=None):
        """Return list of PairedRegions objects w/ mutually conflicting regions

        cm -- ConflictMatrix for this PairedRegions object.

        Mutually conflicting regions form a knot-component as defined
            in Rodland 2006.
        Return value is a list of PairedRegions objects. Each PairedRegions
            object contains mutually conflicting regions (knot-components)
        E.g. region A conflicts with B and B conflicts with A and C, and D
            conflicts with E, and F doesn't conflict with any othe region,
            one group would be A, B and C, the other would be D and E.
            F would not be returned in any group, since it isn't conflicting.
        """
        if cm is None:
            cm = ConflictMatrix(self)
        id_to_pr = self.byId()
        cliques = cm.conflictCliques()
        result = []
        for cl in cliques:
            pr = PairedRegions()
            for i in cl:
                pr.append(id_to_pr[i])
            result.append(pr)
        return result

 
def PairedRegionsFromPairs(pairs):
    """Return PairedRegions object from Pairs

    pairs -- Pairs object, no conflicts allowed

    Result is a list of stretches of perfectly nested base pairs, which means
    [[(m,n), (m+1,n-1), (m+2,n-2),...],[(i,j),(i+1,j-1),...]]
    Base pairs are made directed and sorted before the stretches are
        picked out, so the result will be in order.
    IDs of the regions are set as indices (enumeration of all regions). The
        PairedRegion that starts closest to the 5' end will get ID 0, the next
        ID 1, etc.
    """
    p = Pairs(pairs)
    if not p:
        return PairedRegions()
    if p.hasConflicts():
        raise ValueError("Cannot handle base pair conflicts")
    clean_pairs = p.directed()
    clean_pairs.sort()
    regions = []
    curr_region = []
    pr_id = -1 # paired region ID
    for pair in clean_pairs:
        if not curr_region:
            curr_region.append(pair)
        else:
            x,y = curr_region[-1]
            if pair == (x+1, y-1):
                curr_region.append(pair)
            else:
                pr_id += 1
                regions.append(PairedRegionFromPairs(curr_region, Id=pr_id))
                curr_region = [pair]
    if curr_region: # last block
        pr_id += 1
        regions.append(PairedRegionFromPairs(curr_region, Id=pr_id))
    return PairedRegions(regions)

class ConflictMatrix(object):
    """Stores conflict matrix

    A conflict matrix is a matrix that indicates which PairedRegion objects
        are conflicting. Row and column IDs correspond to Region IDs. If
        two regions are conflicting True is stored, otherwise False is stored.
    """

    def __init__(self, data):
        """Initialize new ConflictMatrix object
        
        data -- either a PairedRegions object or a Pairs object,
            or anything that can be made into a Pairs object
            (e.g. a list of tuples)

        This method sets the Matrix attribute to a Dict2D object containing
            conflict information on the PairedRegions. 
        The input data is either a PairedRegions object or it is made into one.
            A ValueError will be raised when the pairs or regions are
            overlapping. Input data that can't be converted to Pairs
            will lead to downstream errors. Pairs doesn't perform
            any validation.
        Row and column IDs are the Identifiers of the PairedRegion objects.
        The RowOrder and ColumnOrder of the Dict2D are the sorted region IDs.
        """
        if isinstance(data, PairedRegions):
            id_to_pr = data.byId()
        elif isinstance(data, Pairs):
            id_to_pr = PairedRegionsFromPairs(data).byId()
        else: # try to convert to Pairs 
            try:
                d = Pairs(data)
                id_to_pr = PairedRegionsFromPairs(d).byId()
            except:
                raise ValueError("Can't convert data to Pairs")

        # handle the rows and columns in order and set RowOrder and ColOrder
        ro = id_to_pr.keys()
        ro.sort()
        co = id_to_pr.keys()
        co.sort()
        conf = {} # dict of conflicts between blocks
        for id1, bl in id_to_pr.items():
            for id2, bl2 in id_to_pr.items():
                if id2 < id1: # minimize number of calculations
                    continue
                if id1 not in conf:
                    conf[id1] = {}
                if id2 not in conf:
                    conf[id2] = {}
                if id1 == id2:
                    conf[id1][id2] = False
                    conf[id2][id1] = False
                    continue
                is_conflicting = bl.conflicting(bl2)
                conf[id1][id2] = is_conflicting
                conf[id2][id1] = is_conflicting
        self.Matrix = Dict2D(conf, RowOrder=ro, ColOrder=co) # create Dict2D

    def conflictsOf(self, pr_id):
        """Return list of region IDs for regions that conflict with pr_id
        
        pr_id -- row ID in the matrix (ID of paired region)

        Input is ID of a particular region, return value are the IDs of
            all regions that conflict with the given region.
        """
        return [k for k,v in self.Matrix[pr_id].items() if v is True]

    def conflicting(self):
        """Return list of region IDs for conflicting regions
        """
        result = []
        cm = self.Matrix
        for pr_id in cm.RowOrder:
            if contains_true(cm[pr_id].values()):
                result.append(pr_id)
        return result

    def nonConflicting(self):
        """Return list of region IDs for non-conflicting regions
        """
        result = []
        cm = self.Matrix
        for pr_id in cm.RowOrder:
            if not contains_true(cm[pr_id].values()):
                result.append(pr_id)
        return result

    def conflictCliques(self):
        """Return list of lists with IDs of mutually conflicting regions

        See documentation on PairedRegions.conflictCliques for more details.
        """
        cm = self.Matrix

        cliques = []
        seen = {}
        
        for pr_id in cm.RowOrder:
            if pr_id in seen:
                continue
            todo = set([pr_id])
            done = set()
            while todo != done:
                collection = [] # collection of conflicts
                for i in todo:
                    if i in done: # no need to do them twice
                        continue
                    conf = [] # conflicting regions
                    for k,v in cm[i].items():
                        if v is True:
                            conf.append(k)
                    collection.extend(conf) # add conflict to collection
                    done.add(i) # register that i is done
                todo.update(collection) # update todo
            if len(done) > 1:
                cliques.append(list(done))
            for i in done:
                seen[i] = True
        return cliques

# =============================================================================
# SCORING FUNCTIONS FOR DYNAMIC PROGRAMMING APPROACH
# =============================================================================

def num_bps(paired_region):
    """Return number of base pairs (=Length) of paired_region
    
    paired_region -- PairedRegion object
    """
    return paired_region.Length

def hydrogen_bonds(seq):
    """Return function to score a PairedRegion by its hydrogen bonds

    seq -- Sequence object or string

    This method counts the number of hydrogen bonds in Watson-Crick and 
        Wobble base pairs. GC pairs score 3, AU and GU pairs score 2.
    """
    HB_SCORE = {('G','C'): 3, ('C','G'): 3,\
                ('A','U'): 2, ('U','A'): 2,\
                ('G','U'): 2, ('U','G'): 2}

    def apply_to(paired_region):
        """Return score of paired_region by its hydrogen bonds
        
        paired_region -- PairedRegion object

        Scores each base pair in the region by giving each GC base pair 3 
            points and each AU or GU base pair 2 points. Other base pairs
            are ignored and don't add anything to the overall score.
        """
        score = 0
        for up,down in paired_region.Pairs:
            seq_pair = (seq[up],seq[down])
            try:
                score += HB_SCORE[seq_pair]
            except KeyError:
                continue
        return score
    return apply_to

# =============================================================================
# HELPER FUNCTIONS FOR DYNAMIC PROGRAMMING APPROACH
# =============================================================================

def contains_true(i):
    """Return True if input contains True
    
    i -- any object that implements __contains__

    Returns True if True is in the input. Both True and 1 count as True.
    Helper function for ConflictMatrix object.
    """
    try:
        if True in i:
            return True
        return False
    except TypeError: # when i is a string
        return False

def empty_matrix(size):
    """Return square matrix as list of lists of specified size.

    size -- int, number of rows and columns in the matrix.

    This function is a helper function of opt_all.
    Each cell is filled with [PairedRegions()]. This is the initialization
        value needed for a dynamic programming matrix that keeps track
        of all optimal solutions. A solution is a single PairedRegions object.
    """
    if size < 1:
        raise ValueError("The size of the matrix should be at least one")
    result = []
    for i in range(size):
        result.append([])
        for j in range(size):
            result[i].append([PairedRegions()])
    return result

def pick_multi_best(candidates, goal='max'):
    """Return list of unique solutions with a maximum/minimum score.

    candidates -- list of PairedRegions objects

    This function returns a list of all PairedRegions objects that
        have an optimal score (maximum or minimum depending on the goal).
        If the list of candidates is empty, a list containing an
        empty PairedRegions object is returned.
    This function is a helper function of dp_matrix_multi.

    NOTE: PairedRegion IDs must be set. They are checked to avoid
        including unsaturated solutions. Maybe implementation should be
        changed, such that (Start, End, Length) tuples are used as IDs?!
    """
    if not candidates:
        return [PairedRegions()]
    result = []
    best_score = None
    seen = {}
    # Candidates have to be processed in order of length
    can_len = [(c.totalLength(), c) for c in candidates]
    can_len.sort()
    can_len.reverse()
    for l, c in can_len:
        c_ids = tuple(c.sortedIds())
        c_ids_set = set(c_ids)
        if not c or c_ids in seen:
            continue
        this_score = c.totalScore()
        if best_score is None:
            best_score = this_score
            result = [c]
            seen[c_ids] = True
        elif this_score == best_score:
            is_sub = False
            for seen_id in seen:
                if len(c_ids_set) == len(c_ids_set & set(seen_id)):
                    is_sub = True
                    break
            if is_sub:
                seen[c_ids] = True
            else:
                result.append(c)
                seen[c_ids] = True
        elif goal == 'max' and this_score < best_score:
            continue
        elif goal == 'min' and this_score > best_score:
            continue
        else:
            result = [c]
            seen[c_ids] = True
            best_score = this_score

    if not result:
        return [PairedRegions()]
    return result


def dp_matrix_multi(paired_regions, goal='max', scoring_function=num_bps): 
    """Return dynamic programming matrix with top-right half filled

    paired_regions -- PairedRegions object
    goal -- str, 'max' or 'min', if the goal is 'max' the routine
        returns the solutions maximizing the score, if the goal
        is 'min' the solutions with a minimum score are returned.
    scoring_function -- function that can be applied to a PairedRegion
        object and that returns a numerical score. 

    This function fills a matrix that calculates the optimal solution for the 
        pseudoknot-removal problem by storing optimal solutions for smaller
        sub-problems. 

    The number of cells in each DP matrix is the number of given paired
        regions times two, because there is one row and column for each
        start and end point of each region. Only the top-right half of 
        the matrix will be filled. 
    A row index is referred to as i (begin_idx in code), a column index
        is referred to as j (end_idx in code). The matrix is initialized
        on the diagonal (where i==j) with a list containing an empty
        solution (an empty PairedRegions object). A list is used because
        we keep track of all possible optimal choices.
    For each cell (i,j) where j>i we collect all the candidate-solutions
        as follows. 
        ** Add all solutions of the cell to the left, which contains the
        best solutions for the area from start point i to end point j-1.
        ** Add all solutions of the cell to the bottom, which contains
        the best solutions for the area from start point i+1 to end point j.
        ** If start point i and end point j are a start and end point of
        the same region, add all possible solutions from the cell to the
        bottom-left plus this region. The cell to the bottom-left contains
        the optimal solutions for the area from start point i+1 to end
        point j-1.
        ** If the lists of solutions at the cells to the left and bottom
        both contained anything different from the empty solution, we need
        to check two more things:
            ** For each combination of a solution in the left cell and a
            solution in the right cell, calculate the highest end point of
            the solution in the left cell and the lowest start point for the
            bottom cell. In a collection of paired regions, every region has
            an end point, and the highest end point is the largest number in
            the list of all end points. The lowest start value is calculated
            in a similar way.
            ** If the highest end point is lower than the lowest start point,
            it means both solutions are disjoint and can be added to form a
            better solution. Thus, merge the two solutions and add them to the
            list of candidate-solutions.
            ** Otherwise, the solutions are not disjoint, but because of the
            pseudoknots sub-solutions of the two solutions might be combined
            to form a better solution. Create a slider k (splitter in code)
            that runs from the lowest start point minus one to the highest
            end point plus one. For each cell (i,k), (k+1,j) merge all 
            possible solutions and add them to the list of candidate-solutions.
        ** Next, store in cell (i,j) all solutions with an optimal score.
        There might be one solution, or there might be multiple solutions.
        ** Finish the calculation when the top-right cell in the matrix is
        filled. This cell contains the optimal solutions for the given set
        of paired regions.
    """
    if goal not in ['max','min']:
        raise ValueError("goal has to be 'min' or 'max', but is '%s'"%(goal))
    prs = paired_regions
    num_cells = len(prs)*2

    # pre-calculate scores 
    for pr in prs:
        pr.score(scoring_function)

    # create and initialize matrix
    result = empty_matrix(num_cells)
    
    # create some lookup dictionaries
    # enumerated start and end points
    enum_boundaries = prs.enumeratedBoundaries()
    # inverted enumerated start/end points {pr.Start/End: position in list}
    inv_enum_boundaries = prs.invertedEnumeratedBoundaries() 
    pos_to_pr = prs.byStartEnd() # {(pr.Start, pr.End): PairedRegion}
    
    # fill the matrix 
    for end_idx in range(num_cells):
        for begin_idx in range(end_idx-1, -1, -1):
            # look up sequence positions that match indices
            begin_pos = enum_boundaries[begin_idx]
            end_pos = enum_boundaries[end_idx]
            # look up solutions in left and bottom cells 
            left_cell = result[begin_idx][end_idx-1]
            bottom_cell = result[begin_idx+1][end_idx]
            # collect candidates
            candidates = []
            
            # add solutions from the left cell
            for sol in bottom_cell:
                candidates.append(sol)
            # add solutions from the bottom cell
            for sol in left_cell:
                candidates.append(sol)
            # if begin_pos and end_pos are paired:
            # ==> add left_bottom + this region
            if (begin_pos, end_pos) in pos_to_pr:
                this_region = pos_to_pr[(begin_pos, end_pos)]
                bottom_left = result[begin_idx+1][end_idx-1]
                for sol in bottom_left:
                    candidates.append(\
                    PairedRegions(sol+PairedRegions([this_region])))
            
            # if we have a solution in the left and in the bottom cell:
            if left_cell !=[[]] and bottom_cell != [[]]:

                # check whether they can be added or iterate
                for sol1 in left_cell:
                    for sol2 in bottom_cell:
                        he_pos = sol1.highestEnd()
                        he_idx = inv_enum_boundaries[he_pos]
                        ls_pos = sol2.lowestStart()
                        ls_idx = inv_enum_boundaries[ls_pos]
                        
                        # If both solutions are disjoint
                        if he_pos < ls_pos:
                            candidates.append(sol1.merge(sol2))
                        else: # not disjoint ==> iterate
                            for splitter in range(ls_idx-1, he_idx+1):
                                cell_to_left = result[begin_idx][splitter]
                                cell_to_bottom = result[splitter+1][end_idx]
                                if cell_to_bottom == [[]]:
                                    break
                                for sub_sol1 in cell_to_left:
                                    for sub_sol2 in cell_to_bottom:
                                        both = sub_sol1.merge(sub_sol2)
                                        candidates.append(both)
            
            # select all the candidates of maximum length
            best_candidate = pick_multi_best(candidates, goal=goal)
            result[begin_idx][end_idx] = best_candidate

    # return the whole matrix 
    return result

def matrix_solutions(paired_regions, goal='max', scoring_function=num_bps):
    """Return the list of solutions in the top-right cell of the DP matrix

    paired_regions -- PairedRegions object

    This methods fills a dynamic programming matrix (by calling
        dp_matrix_multi) and returns the list of solutions in the 
        top-right cell.
    """
    return dp_matrix_multi(paired_regions, goal=goal,\
        scoring_function=scoring_function)[0][-1]

# =============================================================================
# DYNAMIC PROGRAMMING APPROACH
# =============================================================================

# DP function used to remove pseudoknots
def opt_all(pairs, return_removed=False, goal='max',\
    scoring_function=num_bps):
    """Return a list of pseudoknot-free Pairs objects

    pairs -- Pairs object or list of tuples. One base can only interact
        with one other base, otherwise an error will be raised.
    return_removed -- boolean, if True a list of tuples of
        (nested pairs, removed pairs) will be returned. 
        Default is False --> list of nested pairs only is returned.
    goal -- str, 'max' or 'min', if the goal is 'max' the routine
        returns the solutions maximizing the score, if the goal
        is 'min' the solutions with a minimum score are returned.
    scoring_function -- function that can be applied to a PairedRegion
        object and that returns a numerical score.

    OPTIMIZATION, ALL SOLUTIONS (OA) -- PSEUDOKNOT REMOVAL METHOD

    This method will find all nested structures with an optimal score.
        Since there might be multiple optimal solutions, the retun value
        is always a list. If there is only one solution, the list will
        contain a single element.
    The problem is solved by dynamic programming. For each clique of mutually
        conflicting paired regions a matrix is filled out, and the best
        solutions are added to the result. Non-conflicting regions are part
        of the solution, and don't need to be processed.
    The user can specify the goal (maximize or minimize) and the scoring
        function. If one specifies for example goal='max' and 
        scoring_function=num_bps, the routine finds the nested structures
        with the maximum number of base pairs.
    See documentation of dp_matrix_multi for recursion rules used in the 
        approach.
    """
    if not pairs.hasPseudoknots():
        return [pairs]

    prs = PairedRegionsFromPairs(pairs)
    id_to_bl = prs.byId()
    cm = ConflictMatrix(prs)

    nc_regions = prs.nonConflicting(cm=cm)
    cliques = prs.conflictCliques(cm=cm) 
    
    # basis for all nested structures are the non-conflicting regions
    result = [PairedRegions(nc_regions)]
    # resolve conflicts, store survivors and removed
    for cl in cliques: 
        new_result = []
        best = matrix_solutions(cl, goal=goal,\
            scoring_function=scoring_function)
        for best_sol in best:
            for prev_res in result:
                new_result.append(prev_res.merge(best_sol))
        result = new_result

    if return_removed:
        # collect the removed pairs for each solution
        surviving_ids = []
        for sol in result:
            surviving_ids.append(dict.fromkeys([pr.Id for pr in sol]))
        removed = []
        for sol, surv in zip(result, surviving_ids):
            rem = []
            for pr_id in id_to_bl:
                if pr_id not in surv:
                    rem.extend(id_to_bl[pr_id].Pairs)
            rem.sort()
            removed.append(rem)
        nested = [prs.toPairs() for prs in result]
        return zip(nested, removed)
    
    nested = [prs.toPairs() for prs in result] 
    return nested

# =============================================================================
# MAJORITY OF BASE PAIRS -- CONVENIENCE FUNCTIONS
# =============================================================================

def opt_single_random(pairs, return_removed=False, goal='max',\
    scoring_function=num_bps):
    """Return single pseudoknot-free Pairs object with an optimal score

    pairs -- Pairs object or list of tuples. One base can only interact
        with one other base, otherwise an error will be raised.
    return_removed -- boolean, if True a tuple of (nested pairs, removed pairs)
        will be returned. Default is False --> only nested pairs are returned.
    goal -- str, 'max' or 'min', if the goal is 'max' the routine
        returns the solutions maximizing the score, if the goal
        is 'min' the solutions with a minimum score are returned.
    scoring_function -- function that can be applied to a PairedRegion
        object and that returns a numerical score.

    There might be multiple nested structures with an optimal score.
        This method calculates all of them and returns one at random.
    The user can specify the goal (maximize or minimize) and the scoring
        function. If one specifies for example goal='max' and 
        scoring_function=num_bps, the routine finds the nested structures
        with the maximum number of base pairs.
    """
    nested_structs = opt_all(pairs, return_removed, goal,\
        scoring_function)
    return choice(nested_structs)

def opt_single_property(pairs, return_removed=False, goal='max',\
    scoring_function=num_bps):
    """Return single pseudoknot-free Pairs object with max number of bps

    pairs -- Pairs object or list of tuples. One base can only interact
        with one other base, otherwise an error will be raised.
    return_removed -- boolean, if True a tuple of (nested pairs, removed pairs)
        will be returned. Default is False --> only nested pairs are returned.
    goal -- str, 'max' or 'min', if the goal is 'max' the routine
        returns the solutions maximizing the score, if the goal
        is 'min' the solutions with a minimum score are returned.
    scoring_function -- function that can be applied to a PairedRegion
        object and that returns a numerical score.

    There might be multiple nested structures with an optimal score.
        This method calculates all of them and
        returns the best by examining some properties. The first criterion
        is the number of paired regions in the returned structure, the
        second is the average range of the regions, the third is the 
        average start value of the regions. It returns the structure
        with the minimum value for these three properties. If all properties
        are the same for multiple structures, an error is raises. I believe
        this can't happen, but if it does the behavior can be changed.
    The user can specify the goal (maximize or minimize) and the scoring
        function. If one specifies for example goal='max' and 
        scoring_function=num_bps, the routine finds the nested structures
        with the maximum number of base pairs.
    """
    nested_structs = opt_all(pairs, return_removed, goal,\
        scoring_function)
    lookup = {}
    if return_removed:
        for p, p_rem in nested_structs:
            prs = PairedRegionsFromPairs(p)
            num_regions = len(prs)
            avg_range = average([pr.range() for pr in prs])
            avg_start = average([pr.Start for pr in prs])
            three = (num_regions, avg_range, avg_start)
            if three not in lookup:
                lookup[three] = []
            lookup[three].append((p,p_rem))
    else:
        for p in nested_structs:
            prs = PairedRegionsFromPairs(p)
            num_regions = len(prs)
            avg_range = average([pr.range() for pr in prs])
            avg_start = average([pr.Start for pr in prs])
            three = (num_regions, avg_range, avg_start)
            if three not in lookup:
                lookup[three] = []
            lookup[three].append(p)
    min_key = min(lookup.keys())
    min_value = lookup[min_key]
    if len(min_value) == 1:
        return min_value[0]
    else:
        # believe this can never happen, but just to be sure...
        raise ValueError("Multiple solutions found with equal properties")


# =============================================================================
# CONFLICT-ELIMINATION APPROACHES
# =============================================================================

def find_max_conflicts(conflicting_ids, cm, id_to_pr):
    """Return region ID of the region involved in the most conflicts

    conflicting_ids -- list of PairedRegion IDs
    cm -- ConflictMatrix object
    id_to_pr -- dict of {region ID: PairedRegion}. Result of 
        PairedRegions.byId() method.

    This methods returns the region ID (out of conflicting_ids) involved
        in the most conflicts. If there is a single region with the most
        conflicts, return it. Otherwise compare all regions with the 
        max number of conflicts on their gain. Gain is the length of
        the region minus the cumulative length of all of its conflicting
        regions. Return the one with the minimum gain. If both properties
        are equal, return the region that starts closest to the 3' end.
    """
    number_of_conflicts = {}
    for pr_id in conflicting_ids:
        noc = len(cm.conflictsOf(pr_id))
        if noc not in number_of_conflicts:
            number_of_conflicts[noc] = []
        number_of_conflicts[noc].append(pr_id)
    max_noc = max(number_of_conflicts.keys())
    max_ids = number_of_conflicts[max_noc]
    if len(max_ids) == 1:
        return max_ids[0]
    else:
        len_diffs = {}
        for pr_id in max_ids:
            pr_len = id_to_pr[pr_id].Length
            conf_len = sum([id_to_pr[i].Length for i in cm.conflictsOf(pr_id)])
            diff = pr_len - conf_len
            if diff not in len_diffs:
                len_diffs[diff] = []
            len_diffs[diff].append(pr_id)
        min_ld = min(len_diffs.keys())
        min_ids = len_diffs[min_ld]
        if len(min_ids) == 1:
            return min_ids[0]
        else:
            start_vals = {}
            for pr_id in min_ids:
                start = id_to_pr[pr_id].Start
                start_vals[start] = pr_id
            max_start = max(start_vals.keys())
            return start_vals[max_start]

def find_min_gain(conflicting_ids, cm, id_to_pr):
    """Return region ID of the region with the minimum gain

    conflicting_ids -- list of PairedRegion IDs
    cm -- ConflictMatrix object
    id_to_pr -- dict of {region ID: PairedRegion}. Result of 
        PairedRegions.byId() method.

    This methods returns the region ID (out of conflicting_ids) of the region
        that has the minimum gain. Gain is the length of the region minus
        the cumulative length of all of its conflicting regions. It expresses
        how many base pairs are gained if this region is kept and all
        of its conflicts have to be removed. If its gain is positive, it is 
        favorable to keep this region. If its gain is negative, it is better
        to remove this region and keep its conflicts. If there are multiple
        regions with the minimal gain, the one involved in the most conflicts
        is returned. If both properties are equal, the method returns the
        region that starts closest to the 3' end.
    """
    len_diffs = {}
    for pr_id in conflicting_ids:
        pr_len = id_to_pr[pr_id].Length
        conf_len = sum([id_to_pr[i].Length for i in cm.conflictsOf(pr_id)])
        diff = pr_len - conf_len
        if diff not in len_diffs:
            len_diffs[diff] = []
        len_diffs[diff].append(pr_id)
    min_ld = min(len_diffs.keys())
    min_ids = len_diffs[min_ld]
    if len(min_ids) == 1:
        return min_ids[0]
    else:
        number_of_conflicts = {}
        for pr_id in min_ids:
            noc = len(cm.conflictsOf(pr_id))
            if noc not in number_of_conflicts:
                number_of_conflicts[noc] = []
            number_of_conflicts[noc].append(pr_id)
        max_noc = max(number_of_conflicts.keys())
        max_ids = number_of_conflicts[max_noc]
        if len(max_ids) == 1:
            return max_ids[0]
        else:
            start_vals = {}
            for pr_id in min_ids:
                start = id_to_pr[pr_id].Start
                start_vals[start] = pr_id
            max_start = max(start_vals.keys())
            return start_vals[max_start]

def add_back_non_conflicting(paired_regions, removed):
    """Return new PairedRegions object and new dict of removed regions

    paired_regions -- PairedRegions object
    removed -- dict of {region_id: PairedRegion}

    Helper-function for conflict_elimination. 
    Circular removal might occur in conflict-elimination methods. It means
        that a particular region is removed and later in the process all
        of its conflicts are also removed, which result in an eliminated
        region that doens't conflict with any region in the solution anymore.
    This methods adds removed regions back into the solution if 
        they don't conflict with any region in the solution. The order in 
        which regions are tried to add is from 5' to 3' starting point.
    """
    id_to_pr = paired_regions.byId()
    new_removed = removed.copy()
    
    added = True
    # process removed from 5' to 3'
    order = [(pr.Start, pr.Id) for pr in new_removed.values()]
    order.sort() # from low start value to high start value
    while added:
        added = False
        for start, region_id in order:
            pr1 = new_removed[region_id]
            is_conflicting = False
            for pr2 in id_to_pr.values():
                if pr1.conflicting(pr2):
                    is_conflicting = True
                    new_removed[region_id] = pr1
                    break
            if not is_conflicting:
                id_to_pr[region_id] = pr1
                del new_removed[region_id]
                order = [(pr.Start, pr.Id) for pr in new_removed.values()]
                order.sort() # from low start value to high start value
                added = True
                break

    return PairedRegions(id_to_pr.values()), new_removed


# Conflict-elimination heuristic.
def conflict_elimination(pairs, sel_function, add_back=True,\
    return_removed=False):
    """Return pseudoknot-free Pairs object

    pairs -- Pairs object or list of tuples
    sel_function -- function that takes a list of IDs of conflicting regions,
        a conflict matrix and a dict of {region_id: PairedRegion} and returns
        the ID of a paired region that has to be removed.
    add_back -- boolean, if True regions that are removed but not conflicting
        at the end because of circular removal are added back into the 
        solution. If False, regions are only removed. This choice might result
        in too many regions being removed. Default value is True.
    return_removed -- boolean, if True a tuple of (nested pairs, removed pairs)
        will be returned. Default is False --> only nested pairs are returned.
   
    CONFLICT ELIMINATION -- PSEUDOKNOT REMOVAL METHOD
    EC -- sel_function=find_max_conflicts
    EG -- sel_functino=find_min_gain

    This is the general conflict-elimination function that should
        be used to remove pseudoknots from a knotted RNA structure.
    This algorithm removes paired regions one at the time. The order is 
        in which regions are removed is specified by the selection function.
    Different selection functions can be specified. Selection functions should
        take a list of conflicting IDS, a ConflictMatrix and a dict of
        {Region ID: PairedRegion} as input and they should return a 
        single PairedRegion ID.
    Two selection functions are available: find_max_conflicts and 
        find_min_gain. See their documentation for specifications.
    """
    prs = PairedRegionsFromPairs(pairs)
    id_to_pr = prs.byId()
    cm = ConflictMatrix(prs)
    removed = {}

    conf = cm.conflicting()
    while conf:
        to_remove = sel_function(conf, cm, id_to_pr)
        removed[to_remove] = id_to_pr[to_remove]
        prs.remove(id_to_pr[to_remove])
        id_to_pr = prs.byId()
        cm = ConflictMatrix(prs)
        conf = cm.conflicting()
    # potential circular removal: add regions back in
    if add_back:
        # collect IDs of non-conflicting removed regions
        prs, removed = add_back_non_conflicting(prs, removed)
        
    if return_removed:
        rem = PairedRegions(removed.values()).toPairs()
        return prs.toPairs(), rem
    return prs.toPairs()


# =============================================================================
# INCREMENTAL APPROACHES
# =============================================================================

# Incremental in order (IO) method
def inc_order(pairs, reversed=False, return_removed=False):
    """Return pseudoknot-free Pairs object

    pairs -- Pairs object or list of tuples. One base can only interact
        with one other base, otherwise an error will be raised.
    reversed -- boolean, indication whether the algoritm adds paired
        regions from 5' to 3' start values or from 3' to 5' end values. 
        If False, order is 5' to 3', if True, order is 3' to 5'. Default
        is False.
    return_removed -- boolean, if True a tuple of (nested pairs, removed pairs)
        will be returned. Default is False --> only nested pairs are returned.

    INCREMENTAL IN ORDER (IO) -- PSEUDOKNOT REMOVAL METHOD

    This algorithm treats all the paired regions in order, either starting
    at he 5' end (reversed=F) or at the 3' end (reversed=T). It accepts
    all the non-conflicting regions; If a region conflicts with an already
    added region, it is excluded from the solution.
    """
    prs = PairedRegionsFromPairs(pairs)
    id_to_pr = prs.byId()
    cm = ConflictMatrix(prs)

    if reversed:
        by_pos = [(pr.End, pr) for pr in prs]
        by_pos.sort()
        by_pos.reverse()
    else:
        by_pos = [(pr.Start, pr) for pr in prs]
        by_pos.sort()

    excluded = {}
    result = PairedRegions()
    for pos, pr in by_pos:
        if pr.Id in excluded:
            continue
        result.append(pr)
        for k in cm.conflictsOf(pr.Id):
            excluded[k] = True

    if return_removed:
        removed = Pairs([])
        for pr_id in excluded:
            removed.extend(id_to_pr[pr_id].Pairs)
        removed.sort()
        return result.toPairs(), removed
    return result.toPairs()
    
# Incremental by length (IL) method
def inc_length(pairs, reversed=False, return_removed=False):
    """Return pseudoknot-free Pairs object

    pairs -- Pairs object or list of tuples. One base can only interact
        with one other base, otherwise an error will be raised.
    reversed -- boolean. In case of equal lengths, all paired regions
        are processed from 5' to 3' starting position. 
        If reversed is True, regions are processed from 3' to 5' starting
        position.
    return_removed -- boolean, if True a tuple of (nested pairs, removed pairs)
        will be returned. Default is False --> only nested pairs are returned.

    INCREMENTAL BY LENGTH (IL) -- PSEUDOKNOT REMOVAL METHOD

    This algorithm will process the paired regions from the longest
        to the shortest. In case there are multiple regions of the same length,
        the one on the 5' side is added first if reversed=False (3' side is
        preferred if reversed=True). All paired regions that are conflicting
        with an already-added region are excluded.
    """
    prs = PairedRegionsFromPairs(pairs)
    id_to_pr = prs.byId()
    # create conflict matrix to lookup the conflicts
    cm = ConflictMatrix(prs)

    length_pos_data = {} # dict of {region_length: [(pr.start, pr)]}
    for pr in prs:
        if pr.Length not in length_pos_data:
            length_pos_data[pr.Length] = []
        length_pos_data[pr.Length].append((pr.Start, pr))
    for v in length_pos_data.values():
        v.sort()
        if reversed:
            v.reverse()
    
    excluded = {}
    result = PairedRegions()
   
    lengths = length_pos_data.keys()
    lengths.sort()
    lengths.reverse() # longest regions first

    for pr_len in lengths:
        for pr_start, pr in length_pos_data[pr_len]:
            if pr.Id not in excluded:
                result.append(pr)
                # use the conflict matrix to determine which regions to exclude
                for k in cm.conflictsOf(pr.Id):
                    excluded[k] = True

    if return_removed:
        removed = Pairs([])
        for pr_id in excluded:
            removed.extend(id_to_pr[pr_id].Pairs)
        removed.sort()
        return result.toPairs(), removed
    return result.toPairs()

# Incremental by range (IR) method
def inc_range(pairs, reversed=False, return_removed=False):
    """Return pseudoknot-free Pairs object

    pairs -- Pairs object or list of tuples. One base can only interact
        with one other base, otherwise an error will be raised.
    reversed -- boolean. If reversed is True: in case of two regions
        with the same range the region that starts closest to the 5' side
        is added first. If reversed is False: the region that starts closest
        to the 3' side is added first. Default is False (5' region preferred).
    return_removed -- boolean, if True a tuple of (nested pairs, removed pairs)
        will be returned. Default is False --> only nested pairs are returned.
    
    INCREMENTAL BY RANGE (IR) -- PSEUDOKNOT REMOVAL METHOD

    This algorithm will process the paired regions from the one with the 
        shortest range to the one with the longest range. The range of 
        a region is defined as the distance between the highest upstream
        position and the lowest downstream position (-1); in other words, it
        is the number of unpaired bases in the haripin if this region were
        the only paired region in the structure. 
    In case there are multiple regions of the same length,
        the one that starts closest to the 5' side is added first if 
        reversed=False (starting at the 3' side is preferred if reversed=True).
    All paired regions that are conflicting with an already-added region
        are excluded.
    """
    
    prs = PairedRegionsFromPairs(pairs)
    id_to_pr = prs.byId()
    # create conflict matrix to lookup the conflicts
    cm = ConflictMatrix(prs)

    range_pos_data = {} # dict of {region_range: [(pr.start, pr)]}
    for pr in prs:
        rr = pr.range() # region range
        if rr not in range_pos_data:
            range_pos_data[rr] = []
        range_pos_data[rr].append((pr.Start, pr))
    for v in range_pos_data.values():
        v.sort()
        if reversed:
            v.reverse()

    ranges = range_pos_data.keys()
    ranges.sort()

    result = PairedRegions()
    excluded = {}
    for rr in ranges:
        for pr_start, pr in range_pos_data[rr]:
            if pr.Id not in excluded:
                result.append(pr)
                for k in cm.conflictsOf(pr.Id):
                    excluded[k] = True

    if return_removed:
        removed = Pairs([])
        for pr_id in excluded:
            removed.extend(id_to_pr[pr_id].Pairs)
        removed.sort()
        return result.toPairs(), removed
    return result.toPairs()

# =============================================================================
# NUSSINOV RESTRICTED
# =============================================================================

def nussinov_fill(pairs, size):
    """Return filled dynamic programming search matrix with number of base pairs

    pairs -- Pairs object or list of tuples, should be directed (up,down)
    size -- int, number of rows and columns in the matrix (should be
        at least as much as the highest base paired position)

    Applies Nussinov-Jacobson algorithm restricted to input list of pairs.
    This function records the number of base pairs in the optimal
        (sub)solution.
    """
    bp_dict = dict.fromkeys(pairs)
    m = zeros((size, size), int)
    for j in range(size):
        for i in range(j-1,-1,-1):
            m[i,j] = m[i+1,j] # i unpaired
            if m[i,j-1] > m[i,j]: # j unpaired
                m[i,j] = m[i,j-1]
            if (i,j) in bp_dict and m[i+1,j-1]+1 > m[i,j]: # (i,j) pair
                m[i,j] = m[i+1,j-1]+ 1
            for k in range(i+1,j-1): # bifurcation
                if m[i,k] + m[k+1,j] > m[i,j]:
                    m[i,j] = m[i,k]+m[k+1,j]
    return m

def nussinov_traceback(m, i, j, pairs):
    """Return set of base pairs: nested structure with max number of pairs

    m -- filled DP search matrix
    i -- int, row coordinate where traceback should start, normally 0
    j -- int, column coordinate where traceback should start, normally length-1
    pairs -- Pairs object of list of tuples, should be directed (up, down),
        expect the same list as at the fill stage.

    Traceback procedure of the Nussinov-Jacobson algorithm that returns
        a single solution containing the maximum number of base pairs.
    """
    bp_dict = dict.fromkeys(pairs)
    if m[i,j] == 0: #or if i>=j:
        return set()
    if (i,j) in bp_dict and m[i+1,j-1] + 1 == m[i,j]:
        return set([(i,j)]) | nussinov_traceback(m, i+1, j-1, pairs)
    for k in range(i,j):
        if m[i,j] == m[i,k] + m[k+1,j]:
            return nussinov_traceback(m,i,k,pairs) | \
                nussinov_traceback(m,k+1,j,pairs)

def nussinov_restricted(pairs, return_removed=False):
    """Return nested Pairs object containing the maximum number of pairs

    pairs -- Pairs object or list of tuples [(1,10),(2,9),...]. List
        of base pairs has to be conflict-free (one base can only
        pair with one other base)
    return_removed -- boolean, if True a list of tuples of
        (nested pairs, removed pairs) will be returned. 
        Default is False --> list of nested pairs only is returned.

    This function is a modification of the original Nussinov-Jacobson
        algorithm, restricted to the given list of base pairs.
        It calculated a nested structure containing the maximum
        number of base pairs.

    NOTE: This function is very slow. If have have many base pairs in 
        the list, the size of the matrix grows quickly. We recoomend using
        the opt_all (OA) function for larger problems.
    """
    # pairs have to be conflict-free, directed, and sorted
    p = Pairs(pairs)
    if p.hasConflicts():
        raise ValueError("Cannot handle base pair conflicts")
    p = p.directed()
    p.sort()
    if not p.hasPseudoknots():
        nested = p  # if not Pseudoknots: return structure
    else:
        paired_positions = [x for x,y in p] + [y for x,y in p]
        paired_positions.sort()
        if max(paired_positions) > 200:
            # if the 'sequence' is longer than 200, map the numbers
            # to remove the 'unpaired' positions. Avoid waste of space
            mapped_back = dict([(x,y) for x,y in enumerate(paired_positions)])
            mapped = dict([(y,x) for x,y in enumerate(paired_positions)])
            mapped_pairs = []
            for x,y in pairs:
                mapped_pairs.append((mapped[x],mapped[y]))
            max_idx = len(paired_positions)+1
            m = nussinov_fill(mapped_pairs, size=max_idx)
            t = nussinov_traceback(m, 0, max_idx-1, mapped_pairs)
            nested = []
            for x,y in t:
                nested.append((mapped_back[x], mapped_back[y]))
        else:
            # sequence short enough, fill matrix directly
            max_idx = max(filter(None, paired_positions))+1
            m = nussinov_fill(p, size=max_idx)
            nested = nussinov_traceback(m, 0, max_idx-1, p)
            nested = Pairs(list(nested))
        nested.sort()
    if return_removed:
        removed = Pairs([])
        for bp in p:
            if bp not in nested:
                removed.append(bp)
        return nested, removed
    else:
        return nested



if __name__ == "__main__":
    pass