+ changed check for intial setting of old_w. (matrix != None does not work in cvxopt)
[qpalma.git] / qpalma / Plif.py
1 # This program is free software; you can redistribute it and/or modify
2 # it under the terms of the GNU General Public License as published by
3 # the Free Software Foundation; either version 2 of the License, or
4 # (at your option) any later version.
5 #
6 # Written (W) 2008 Fabio De Bona
7 # Copyright (C) 2008 Max-Planck-Society
8
9 import math
10 import pdb
11
12 from numpy.matlib import zeros,isinf
13
14 import QPalmaDP
15
16
17 class Plf:
18 """
19 means piecewise linear function
20 """
21
22 def __init__(self,num=None):
23 if num == None:
24 self.len = 0
25 self.limits = []
26 self.penalties = []
27 self.transform = ''
28 self.name = ''
29 self.max_len = 0
30 self.min_len = 0
31 else:
32 self.len = num
33 self.limits = linspace(0,40,num)
34 self.penalties = [0]*num
35 self.transform = ''
36 self.name = ''
37 self.max_len = 0
38 self.min_len = 0
39
40
41 def convert2SWIG(self):
42 ps = QPalmaDP.penalty_struct()
43
44 ps.len = len(self.limits)
45
46 ps.limits = QPalmaDP.createDoubleArrayFromList(self.limits)
47 ps.penalties = QPalmaDP.createDoubleArrayFromList(self.penalties)
48
49 ps.max_len = self.max_len
50 ps.min_len = self.min_len
51 ps.transform = 0
52 ps.name = self.name
53
54 return ps
55
56
57 def base_coord(b1,b2):
58 b1 = b1.lower()
59 b2 = b2.lower()
60 assert b1 in ['a','c','g','t','n']
61 assert b2 in ['a','c','g','t','n']
62 cb = {'a':0, 'c':1, 'g':2, 't':3, 'n':4}
63
64 b1 = cb[b1]
65 b2 = cb[b2]
66
67 return b1*5+b2
68
69 def linspace(a,b,n):
70 intervalLength = b-a
71 stepSize = 1.0*intervalLength / (n-1)
72
73 interval = [0]*n
74 interval[0] = a
75 interval[-1] = b
76 for i in range(1,n-1):
77 interval[i] = a+(i*stepSize)
78
79 return interval
80
81 def logspace(a,b,n):
82 interval = [0]*n
83 begin = 10.0**a
84 end = 10.0**b
85 intervalSize = 1.0*(b-a)/(n-1)
86 interval[0] = begin
87 interval[-1] = end
88
89 for i in range(1,n-1):
90 interval[i] = 10**(a+i*intervalSize)
91
92 return interval
93
94 def log10(x):
95 return math.log(x)/math.log(10)
96
97
98 def penalty_lookup_new(penalty_struct, value):
99
100 limits = penalty_struct.limits
101 penalties = penalty_struct.penalties
102
103 if penalty_struct.transform == 'log':
104 value = math.log(value)
105
106 elif penalty_struct.transform == 'log(+3)':
107 value = math.log(value+3)
108
109 elif penalty_struct.transform == 'log(+1)':
110 value = math.log(value+1)
111
112 elif penalty_struct.transform == '':
113 pass
114
115 elif penalty_struct.transform == '(+3)':
116 value += 3
117
118 elif penalty_struct.transform == 'mod3':
119 value = value%3
120
121 else:
122 assert False,'unknown penalty transform'
123
124 count = len([elem for elem in limits if elem <= value])
125
126 if count == 0:
127 pen = penalties[0]
128 elif count == len(limits):
129 pen=penalties[-1]
130 else:
131 pen = (penalties[count]*(value-limits[count-1]) + penalties[count-1]\
132 *(limits[count]-value)) / (limits[count]-limits[count-1])
133
134 return pen
135
136
137 def compute_donacc(donor_supp, acceptor_supp, d, a):
138
139 assert(len(donor_supp)==len(acceptor_supp))
140 size = len(donor_supp)
141
142 donor = [0.0]*size
143 acceptor = [0.0]*size
144
145 for idx in range(size):
146 if isinf(donor_supp[idx]):
147 donor[idx] = donor_supp[idx]
148 else:
149 donor[idx] = penalty_lookup_new(d, donor_supp[idx])
150
151 if isinf(acceptor_supp[idx]):
152 acceptor[idx] = acceptor_supp[idx]
153 else:
154 acceptor[idx] = penalty_lookup_new(a,acceptor_supp[idx])
155
156 return donor,acceptor
157
158