#! /usr/bin/env python

# numLocOpt.py
# P.Clote

# Computes the number of locally optimal sec str, assuming
# every base can base pair with every other base. Used for
# combinatorics study of locally optimal secondary structures.

import sys,math
THRESHOLD = 1
R         =  0.0019872370936902486 
T         = 310

def maxNumBasePairs(n):
  return (n-THRESHOLD)/2


def locOpt(n):  
  #D[(m,k)] is number of loc opt str on m having k base pairs
  #M(n,k) num of loc opt on m having k base pairs with NO visible positions
  D = {}; M = {}  
  for m in range(THRESHOLD+2):
    M[(m,0)] = 0L; D[(m,0)] = 1L
    for k in range(1,maxNumBasePairs(m)+1): 
      M[(m,k)] = 0L; D[(m,k)] = 0L
  M[(0,0)] = 1L
  for m in range(THRESHOLD+2,n+1):
    #COMPUTE M
    M[(m,0)] = 0L  #once m>threshold+1, there are no loc opt with 0 bp
    for k in range(1,maxNumBasePairs(m)+1): 
      sum = 0L
      if k-1 <= maxNumBasePairs(m-2):
        sum += D[(m-2,k-1)]  #Case 1: (1,m) is base pair
      for i in range(2,m-THRESHOLD): #Case 3: (i,m) is base pair 
        for k0 in range(0,maxNumBasePairs(i-1)+1):
          #k0 many base pairs in sequence 1,...,i-1
          k1 = k-k0-1  #want k0+k1+1 = k many base pairs
          if k1<0 or k1>maxNumBasePairs(m-i-1): continue
          sum += M[(i-1,k0)]*D[(m-i-1,k1)]
      M[(m,k)] = sum
    #COMPUTE D
    D[(m,0)] = 0L  #once m>threshold+1, there are no loc opt with 0 bp
    for k in range(1,maxNumBasePairs(m)+1): 
      sum = 0L
      if k-1 <= maxNumBasePairs(m-2):
        sum += D[(m-2,k-1)]  #Case 1: (1,m) is base pair
      for r in range(m-THRESHOLD-1,m):  #Case 2: m not base paired
        #Case 2: r+1,...,m unpaired, for some m > r >= m-THRESHOLD-1
        #Note that if r<m-THRESHOLD-1 and r+1,...,m unpaired then S not loc opt
        if k>maxNumBasePairs(r):  continue
        sum += M[(r,k)]  
      for i in range(2,m-THRESHOLD): #Case 3: (i,m) is base pair 
        for k0 in range(0,maxNumBasePairs(i-1)+1):
          #k0 many base pairs in sequence 1,...,i-1
          k1 = k-k0-1  #want k0+k1+1 = k many base pairs
          if k1<0 or k1>maxNumBasePairs(m-i-1): continue
          sum += D[(i-1,k0)]*D[(m-i-1,k1)]
      D[(m,k)] = sum
  return D

def partFunLocOpt(n,D):  
  E = {}; Z = 0L
  keys = []
  for k in range(maxNumBasePairs(n)+1):
    if (n,k) in D.keys(): keys.append(k)
  for k in keys: Z += D[(n,k)] * math.exp(-k/(R*T))
  for k in keys: E[(n,k)] = D[(n,k)] * math.exp(-k/(R*T)) / Z
  return E

def main(n):
  D = locOpt(n)
  E = partFunLocOpt(n,D)
  print "\nk\tNum loc optimal with k basepairs"
  for k in range(maxNumBasePairs(n)+1):
    if D[(n,k)]>0: print "%d\t%s" % (k,D[(n,k)])
  print "\nk\tNum k-loc optimal"
  for k in range(maxNumBasePairs(n),-1,-1):
    k0 = maxNumBasePairs(n)-k
    if D[(n,k)]>0: print "%d\t%s" % (k0,D[(n,k)])
  print "\nk\tPartition function for num loc optimal with k basepairs"
  for k in range(maxNumBasePairs(n)+1):
    if E[(n,k)]>0: print "%d\t%s" % (k,E[(n,k)])
  print "\nk\tPartition function for num k-loc optimal"
  for k in range(maxNumBasePairs(n),-1,-1):
    k0 = maxNumBasePairs(n)-k
    if E[(n,k)]>0: print "%d\t%s" % (k0,E[(n,k)])
  print "\nk\tEnsemble free energy for num k-loc optimal"
  for k in range(maxNumBasePairs(n),-1,-1):
    k0 = maxNumBasePairs(n)-k
    if E[(n,k)]>0: print "%d\t%f" % (k0,math.log(float(E[(n,k)])))
  """
  """
  
  
  

if __name__ == '__main__':
  if len(sys.argv)!= 2:
    print "Usage: %s integer" % sys.argv[0]
    sys.exit(1)
  main(int(sys.argv[1]))