-
Notifications
You must be signed in to change notification settings - Fork 0
/
constrained_softmax.py
119 lines (107 loc) · 3.05 KB
/
constrained_softmax.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
import numpy as np
import pdb
def constrained_softmax(z, u):
z -= np.mean(z)
q = np.exp(z)
active = np.ones(len(u))
mass = 0.
p = np.zeros(len(z))
while True:
inds = active.nonzero()[0]
p[inds] = q[inds] * (1. - mass) / sum(q[inds])
found = False
#import pdb; pdb.set_trace()
for i in inds:
if p[i] > u[i]:
p[i] = u[i]
mass += u[i]
found = True
active[i] = 0
if not found:
break
#print mass
#print active
return p, active, mass
def constrained_softmax_matrix(Z, U):
raise NotImplementedError
def gradient_constrained_softmax_old(z, u, dp, p, active, mass):
n = len(z)
inds = active.nonzero()[0]
Jz = np.zeros((n, n))
Ju = np.zeros((n, n))
for i in xrange(n):
for j in xrange(n):
if active[i]:
if active[j]:
Jz[i, j] = -p[i] * p[j] / (1. - mass)
else:
Ju[i, j] = -p[i] / (1. - mass)
if active[i]:
Jz[i, i] += p[i]
else:
Ju[i, i] = 1.
#print 'Jz: ', Jz
#print 'Ju: ', Ju
dz = Jz.transpose().dot(dp)
du = Ju.transpose().dot(dp)
#import pdb; pdb.set_trace()
return dz, du
def gradient_constrained_softmax(z, u, dp, p, active, mass):
n = len(z)
dp_av = sum(active * p * dp) / (1. - mass)
# inds = active.nonzero()[0]
# dp_av = p[inds].dot(dp[inds]) / (1. - mass)
dz = active * p * (dp - dp_av)
du = (1 - active) * (dp - dp_av)
return dz, du
def numeric_gradient_constrained_softmax(z, u, dp, p, active, mass):
epsilon = 1e-6
n = len(z)
Jz = np.zeros((n, n))
Ju = np.zeros((n, n))
for j in xrange(n):
z1 = z.copy()
z2 = z.copy()
z1[j] -= epsilon
z2[j] += epsilon
p1, _, _ = constrained_softmax(z1, u)
p2, _, _ = constrained_softmax(z2, u)
Jz[:, j] = (p2 - p1) / (2*epsilon)
#import pdb; pdb.set_trace()
u1 = u.copy()
u2 = u.copy()
u1[j] -= epsilon
u2[j] += epsilon
p1, _, _ = constrained_softmax(z, u1)
p2, _, _ = constrained_softmax(z, u2)
Ju[:, j] = (p2 - p1) / (2*epsilon)
#print 'Jz_: ', Jz
#print 'Ju_: ', Ju
dz = Jz.transpose().dot(dp)
du = Ju.transpose().dot(dp)
#import pdb; pdb.set_trace()
return dz, du
if __name__ == "__main__":
n = 6
z = np.random.randn(n)
uc = 0.5*np.random.rand(n)
pdb.set_trace()
#z = np.array([1.,1.,1.,1.,1.,1.])
#uc = np.array([.15,.167,1.,1.,1.,1.])
print sum(uc)
print z
print uc
p, active, mass = constrained_softmax(z, uc)
print p
print sum(p)
dp = np.random.randn(len(z))
pdb.set_trace()
dz, du = gradient_constrained_softmax(z, uc, dp, p, active, mass)
print dp
print dz
print du
dz_, du_ = numeric_gradient_constrained_softmax(z, uc, dp, p, active, mass)
print dz_
print du_
print np.linalg.norm(dz - dz_)
print np.linalg.norm(du - du_)