#! /usr/bin/env python

# expNumHairpinsInSaturated.py
# P.Clote

# Program computes expected number of hairpins in saturated structures, 
# according to recurrence relations in paper submitted by
# Clote, Kranakis and Krizanc to Bulletin of Math Biol in Aug 2013.

import sys,math,random

THRESHOLD = 1

def asymptoticValue(n,txt):
  if txt=="sat":
    return 0.123194*n
  elif txt=="all":
    return 0.105573 * n

def numSatStrHairpin(n):
  #H[m][k] is number of sat str on m having k hairpins
  #J[m][k] is number of sat str on m having k hairpins, with no visible bases
  H = {}; J = {}
  for m in range(n+1):
    H[m]={}; J[m]={}
    for k in range(n+1):
      H[m][k] = 0L; J[m][k] = 0L
  H[0][0] = 1L; H[1][0] = 1L; H[2][0] = 1L; H[3][1] = 1L;
  J[0][0] = 1L; J[3][1] = 1L
  for m in range(THRESHOLD+3,n+1): 
    for k in range(1,n+1): 
      sumH  = 0L;   sumJ = 0L #sumH will be H[m][k], sumJ will be J[m][k]
      sumH += (J[m-1][k] + J[m-2][k])
      for r in range(1,m-THRESHOLD): #r in [1,...,m-2]
        #sat str having k-1 hairpins in left region plus bp (r,k)
        sumH += H[r-1][k-1]*H[m-r-1][0]
        sumJ += J[r-1][k-1]*H[m-r-1][0]
        for s in range(k): #s in [0,...,k-1]
          #sat str having s hairpins in left and k-s in right regions
          #print "m,k,r,s,r-1,s,m-r-1,k-s"
          #print m,k,r,s,r-1,s,m-r-1,k-s
          sumH += H[r-1][s]*H[m-r-1][k-s]
          sumJ += J[r-1][s]*H[m-r-1][k-s]
      H[m][k] = sumH
      J[m][k] = sumJ
  return H,J

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 checkSums(n):
  #check that we correctly compute number of sat str and sat str with no vis.
  A,B = saturatedSecStrRecRel(n)
  R,S = saturatedSecStrGrammar(n)
  H,J = numSatStrHairpin(n)
  for m in range(1,n+1):
    avgH = 0.0; avgJ = 0.0
    for k in range(0,n+1):
      avgH += H[m][k]; avgJ += J[m][k]
    print "%d\t%d\t%d\t%d\t%d\t%d\t%d" % (m,avgH,A[m],R[m],avgJ,B[m],S[m])



def main(n):
  A,B = saturatedSecStrRecRel(n)
#  R,S = saturatedSecStrGrammar(n)
  H,J = numSatStrHairpin(n)
  avgHdict = {}; avgJdict = {}
  for m in range(1,n+1):
    avgH = 0.0; avgJ = 0.0
    for k in range(1,n+1):
      avgH += k*H[m][k]; avgJ += k*J[m][k]
    if A[m]==0:
      avgH = 0
    else:
      avgH = avgH/A[m]
    if B[m]==0:
      avgJ = 0
    else:
      avgJ = avgJ/B[m]  
    avgHdict[m]=avgH; avgJdict[m]=avgJ
    val = asymptoticValue(m,"sat")
    print "%d\t%.8f\t%.8f" % (m,avgH/m,val/m)
  return avgHdict,avgJdict  

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