#! /usr/bin/env python # numSaturatedJCB.py # P.Clote # Program computes number of saturated structures, according # to recurrence relations on page 1646-47 of P. Clote # "Combinatorics of saturated secondary structures of RNA", # J Comp Biol. 13(9):1640-1657 (2006) import sys,math,random THRESHOLD = 1 def asymptoticValue(n): # return 1.07427068741 * n**(-3/2) * 0.424687310420272**(-n) return 1.07427068741 * n**(-3/2) * 2.35467360447**n def saturatedSecStrRecRel(n): #A[m] is number of sat str on m #B[m] is number of sat str on m with no visible bases A = {}; B = {} for m in range(THRESHOLD+2): A[m] = 1L; B[m] = 0L B[0] = 1L for m in range(THRESHOLD+2,n+1): sumA = 0L; sumB = 0L for k in range(1,m-THRESHOLD): sumA += A[k-1]*A[m-k-1] sumB += B[k-1]*A[m-k-1] B[m] = sumB for k in range(1,THRESHOLD+2): sumA += B[m-k] A[m] = sumA B[m] = sumB return A,B def saturatedSecStrGrammar(n): #S[m] is number of non-empty sat str on m #R[m] is number of non-empty sat str on m with no visible points R = {}; S = {} S[1] = 1; S[2]=1; R[1] = 0; R[2] = 0; for m in range(1,THRESHOLD+2): R[m] = 0L S[m] = 1L #empty structure R[THRESHOLD+2] = 1L #(*) S[THRESHOLD+2] = 1L #(*) but not empty structure for m in range(THRESHOLD+3,n+1): S[m] = R[m-1] + R[m-2] + S[m-2] for k in range(1,m-THRESHOLD-1): #k in [1,m-THRESHOLD-2] S[m] += S[k]*S[m-k-2] R[m] = S[m-2] for k in range(1,m-THRESHOLD-1): #k in [1,m-THRESHOLD-2] R[m] += R[k]*S[m-k-2] return R,S def main(n): A,B = saturatedSecStrRecRel(n) R,S = saturatedSecStrGrammar(n) for m in range(1,n+1): val = asymptoticValue(m) # print "%d\t%d\t%d" % (m,S[m],R[m]) print "%d\t%d\t%.f\t%.2f\t%d\t%.2f" % (m,S[m],val,S[m]/val,A[m],A[m]/val) if __name__ == '__main__': if len(sys.argv)!= 2: print "Usage: %s integer" % sys.argv[0] sys.exit(1) main(int(sys.argv[1]))