[MRG] Adding greenkhorn #66
Conversation
arakotom
commented
Sep 24, 2018
arakotom
commented
- added greenkhorn algorithm in bregman.py
- added novel option for resolution
- modified unit test in test_variants
ot/bregman.py
Outdated
| >>> a=[.5,.5] | ||
| >>> b=[.5,.5] | ||
| >>> M=[[0.,1.],[1.,0.]] | ||
| >>> ot.sinkhorn(a,b,M,1) |
ot/bregman.py
Outdated
| m = b.shape[0] | ||
|
|
||
| # Next 3 lines equivalent to K= np.exp(-M/reg), but faster to compute | ||
| K = np.empty(M.shape, dtype=M.dtype) |
ot/bregman.py
Outdated
| np.divide(M, -reg, out=K) | ||
| np.exp(K, out=K) | ||
|
|
||
| u = np.ones(n)/n |
ot/bregman.py
Outdated
| np.exp(K, out=K) | ||
|
|
||
| u = np.ones(n)/n | ||
| v = np.ones(m)/m |
ot/bregman.py
Outdated
|
|
||
| u = np.ones(n)/n | ||
| v = np.ones(m)/m | ||
| G = np.diag(u)@K@np.diag(v) |
There was a problem hiding this comment.
use broadcasting to avoid filling diagonal matrices
G = u[:, np.newaxis] * K * v[, np.newaxis]
ot/bregman.py
Outdated
| G[:,i_2] = u*K[:,i_2]*v[i_2] | ||
| #aviol = (G@one_m - a) | ||
| #aviol_2 = (G.T@one_n - b) | ||
| viol = viol + ( -old_v + v[i_2])*K[:,i_2]*u |
test/test_bregman.py
Outdated
|
|
||
| # check values | ||
| np.testing.assert_allclose(G0, Gs, atol=1e-05) | ||
| np.testing.assert_allclose(G0, Ges, atol=1e-05) | ||
| np.testing.assert_allclose(G0, Gerr) | ||
|
|
||
| np.testing.assert_allclose(G0, G_green, atol = 1e-32) |
There was a problem hiding this comment.
a pep8 checker would tell you but you should not put spaces around = in function signatures. It's to visually distinguish what is a function parameter from a variable assignment.
test/test_bregman.py
Outdated
|
|
||
| # check values | ||
| np.testing.assert_allclose(G0, Gs, atol=1e-05) | ||
| np.testing.assert_allclose(G0, Ges, atol=1e-05) | ||
| np.testing.assert_allclose(G0, Gerr) | ||
|
|
||
| np.testing.assert_allclose(G0, G_green, atol = 1e-32) | ||
| print(G0,G_green) |
There was a problem hiding this comment.
and you should always put a space after a ,
ot/bregman.py
Outdated
| one_n = np.ones(n) | ||
| one_m = np.ones(m) | ||
| viol = G@one_m - a | ||
| viol_2 = G.T@one_n - b |
There was a problem hiding this comment.
here to allocate arrays of ones to compute sum of rows and columns. I would just use np.sum(..., axis=)
ot/bregman.py
Outdated
| log['u'] = u | ||
| log['v'] = v | ||
|
|
||
| while i < numItermax and stopThr_val > stopThr: |
There was a problem hiding this comment.
rather than using a while you could use a for loop. For optim solvers I tend to do:
for i in range(numItermax):
...
if stopping condition satisfied do:
break
else:
print("Solver did not converge")
so you can easily print a message when you did not converge.
ot/bregman.py
Outdated
|
|
||
| """ | ||
|
|
||
| i = 0 |
There was a problem hiding this comment.
indeed, thougt I got rid of all of them...
