+ moved penalty lookup and compute_doncc to Plif module
[qpalma.git] / qpalma / Plif.py
1 #!/usr/bin/env python
2 # -*- coding: utf-8 -*-
3
4 import math
5 import pdb
6
7 from numpy.matlib import zeros,isinf
8
9 import QPalmaDP
10
11
12 class Plf:
13 """
14 means piecewise linear function
15 """
16
17 def __init__(self,num=None):
18 if num == None:
19 self.len = 0
20 self.limits = []
21 self.penalties = []
22 self.transform = ''
23 self.name = ''
24 self.max_len = 0
25 self.min_len = 0
26 else:
27 self.len = num
28 self.limits = linspace(0,40,num)
29 self.penalties = [0]*num
30 self.transform = ''
31 self.name = ''
32 self.max_len = 0
33 self.min_len = 0
34
35
36 def convert2SWIG(self):
37 ps = QPalmaDP.penalty_struct()
38
39 ps.len = len(self.limits)
40
41 ps.limits = QPalmaDP.createDoubleArrayFromList(self.limits)
42 ps.penalties = QPalmaDP.createDoubleArrayFromList(self.penalties)
43
44 ps.max_len = self.max_len
45 ps.min_len = self.min_len
46 ps.transform = 0
47 ps.name = self.name
48
49 return ps
50
51
52 def base_coord(b1,b2):
53 b1 = b1.lower()
54 b2 = b2.lower()
55 assert b1 in ['a','c','g','t','n']
56 assert b2 in ['a','c','g','t','n']
57 cb = {'a':0, 'c':1, 'g':2, 't':3, 'n':4}
58
59 b1 = cb[b1]
60 b2 = cb[b2]
61
62 return b1*5+b2
63
64 def linspace(a,b,n):
65 intervalLength = b-a
66 stepSize = 1.0*intervalLength / (n-1)
67
68 interval = [0]*n
69 interval[0] = a
70 interval[-1] = b
71 for i in range(1,n-1):
72 interval[i] = a+(i*stepSize)
73
74 return interval
75
76 def logspace(a,b,n):
77 interval = [0]*n
78 begin = 10.0**a
79 end = 10.0**b
80 intervalSize = 1.0*(b-a)/(n-1)
81 interval[0] = begin
82 interval[-1] = end
83
84 for i in range(1,n-1):
85 interval[i] = 10**(a+i*intervalSize)
86
87 return interval
88
89 def log10(x):
90 return math.log(x)/math.log(10)
91
92
93 def penalty_lookup_new(penalty_struct, value):
94
95 limits = penalty_struct.limits
96 penalties = penalty_struct.penalties
97
98 if penalty_struct.transform == 'log':
99 value = math.log(value)
100
101 elif penalty_struct.transform == 'log(+3)':
102 value = math.log(value+3)
103
104 elif penalty_struct.transform == 'log(+1)':
105 value = math.log(value+1)
106
107 elif penalty_struct.transform == '':
108 pass
109
110 elif penalty_struct.transform == '(+3)':
111 value += 3
112
113 elif penalty_struct.transform == 'mod3':
114 value = value%3
115
116 else:
117 assert False,'unknown penalty transform'
118
119 count = len([elem for elem in limits if elem <= value])
120
121 if count == 0:
122 pen = penalties[0]
123 elif count == len(limits):
124 pen=penalties[-1]
125 else:
126 pen = (penalties[count]*(value-limits[count-1]) + penalties[count-1]\
127 *(limits[count]-value)) / (limits[count]-limits[count-1])
128
129 return pen
130
131
132 def compute_donacc(donor_supp, acceptor_supp, d, a):
133
134 assert(len(donor_supp)==len(acceptor_supp))
135 size = len(donor_supp)
136
137 donor = [0.0]*size
138 acceptor = [0.0]*size
139
140 for idx in range(size):
141 if isinf(donor_supp[idx]):
142 donor[idx] = donor_supp[idx]
143 else:
144 donor[idx] = penalty_lookup_new(d, donor_supp[idx])
145
146 if isinf(acceptor_supp[idx]):
147 acceptor[idx] = acceptor_supp[idx]
148 else:
149 acceptor[idx] = penalty_lookup_new(a,acceptor_supp[idx])
150
151 return donor,acceptor
152
153