-
-
Notifications
You must be signed in to change notification settings - Fork 2k
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Fix sparse scaling matrix in quadpotential
- Loading branch information
Showing
2 changed files
with
184 additions
and
24 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,152 @@ | ||
import functools | ||
import numpy as np | ||
import scipy.sparse | ||
import theano.tensor as tt | ||
import theano | ||
|
||
from pymc3.step_methods.hmc import quadpotential | ||
|
||
from nose.tools import raises, assert_raises | ||
from nose.plugins.skip import SkipTest | ||
|
||
if quadpotential.chol_available: | ||
from scipy import sparse | ||
from sksparse import cholmod | ||
|
||
|
||
def require_sparse(f): | ||
@functools.wraps(f) | ||
def inner(): | ||
if not quadpotential.chol_available: | ||
raise SkipTest("Test requires sksparse.cholmod") | ||
f() | ||
return inner | ||
|
||
|
||
@raises(quadpotential.PositiveDefiniteError) | ||
def test_elemwise_posdef(): | ||
scaling = np.array([0, 2, 3]) | ||
quadpotential.quad_potential(scaling, True, True) | ||
|
||
|
||
@raises(quadpotential.PositiveDefiniteError) | ||
def test_elemwise_posdef2(): | ||
scaling = np.array([0, 2, 3]) | ||
quadpotential.quad_potential(scaling, True, False) | ||
|
||
|
||
def test_elemwise_velocity(): | ||
scaling = np.array([1, 2, 3]) | ||
x_ = np.ones_like(scaling) | ||
x = tt.vector() | ||
x.tag.test_value = x_ | ||
pot = quadpotential.quad_potential(scaling, True, False) | ||
v = theano.function([x], pot.velocity(x)) | ||
assert np.allclose(v(x_), scaling) | ||
pot = quadpotential.quad_potential(scaling, True, True) | ||
v = theano.function([x], pot.velocity(x)) | ||
assert np.allclose(v(x_), 1/scaling) | ||
|
||
|
||
def test_elemwise_energy(): | ||
scaling = np.array([1, 2, 3]) | ||
x_ = np.ones_like(scaling) | ||
x = tt.vector() | ||
x.tag.test_value = x_ | ||
pot = quadpotential.quad_potential(scaling, True, False) | ||
energy = theano.function([x], pot.energy(x)) | ||
assert np.allclose(energy(x_), 0.5 * scaling.sum()) | ||
pot = quadpotential.quad_potential(scaling, True, True) | ||
energy = theano.function([x], pot.energy(x)) | ||
assert np.allclose(energy(x_), 0.5 * (1/scaling).sum()) | ||
|
||
|
||
def test_equal_diag(): | ||
np.random.seed(42) | ||
for _ in range(3): | ||
diag = np.random.rand(5) | ||
x_ = np.random.randn(5) | ||
x = tt.vector() | ||
x.tag.test_value = x_ | ||
pots = [ | ||
quadpotential.quad_potential(diag, False, False), | ||
quadpotential.quad_potential(1. / diag, True, False), | ||
quadpotential.quad_potential(np.diag(diag), False, False), | ||
quadpotential.quad_potential(np.diag(1. / diag), True, False), | ||
] | ||
if quadpotential.chol_available: | ||
diag_ = sparse.csc_matrix(np.diag(1. / diag)) | ||
pots.append(quadpotential.quad_potential(diag_, True, False)) | ||
|
||
v = np.diag(1. / diag).dot(x_) | ||
e = x_.dot(np.diag(1. / diag).dot(x_)) / 2 | ||
for pot in pots: | ||
v_function = theano.function([x], pot.velocity(x)) | ||
e_function = theano.function([x], pot.energy(x)) | ||
assert np.allclose(v_function(x_), v) | ||
assert np.allclose(e_function(x_), e) | ||
|
||
|
||
def test_equal_dense(): | ||
np.random.seed(42) | ||
for _ in range(3): | ||
cov = np.random.rand(5, 5) | ||
cov += cov.T | ||
cov += 10 * np.eye(5) | ||
inv = np.linalg.inv(cov) | ||
assert np.allclose(inv.dot(cov), np.eye(5)) | ||
x_ = np.random.randn(5) | ||
x = tt.vector() | ||
x.tag.test_value = x_ | ||
pots = [ | ||
quadpotential.quad_potential(cov, False, False), | ||
quadpotential.quad_potential(inv, True, False), | ||
] | ||
if quadpotential.chol_available: | ||
pots.append(quadpotential.quad_potential(cov, False, False)) | ||
|
||
v = np.linalg.solve(cov, x_) | ||
e = x_.dot(v) / 2 | ||
for pot in pots: | ||
v_function = theano.function([x], pot.velocity(x)) | ||
e_function = theano.function([x], pot.energy(x)) | ||
assert np.allclose(v_function(x_), v) | ||
assert np.allclose(e_function(x_), e) | ||
|
||
|
||
def test_random_diag(): | ||
d = np.arange(10) + 1 | ||
pots = [ | ||
quadpotential.quad_potential(d, True, False), | ||
quadpotential.quad_potential(1./d, False, False), | ||
quadpotential.quad_potential(np.diag(d), True, False), | ||
quadpotential.quad_potential(np.diag(1./d), False, False), | ||
] | ||
if quadpotential.chol_available: | ||
d_ = scipy.sparse.csc_matrix(np.diag(d)) | ||
pot = quadpotential.quad_potential(d, True, False) | ||
pots.append(pot) | ||
for pot in pots: | ||
vals = np.array([pot.random() for _ in range(10000)]) | ||
assert np.allclose(vals.std(0), np.sqrt(1./d), atol=0.1) | ||
|
||
|
||
def test_random_dense(): | ||
np.random.seed(42) | ||
for _ in range(3): | ||
cov = np.random.rand(5, 5) | ||
cov += cov.T | ||
cov += 10 * np.eye(5) | ||
inv = np.linalg.inv(cov) | ||
assert np.allclose(inv.dot(cov), np.eye(5)) | ||
|
||
pots = [ | ||
quadpotential.QuadPotential(cov), | ||
quadpotential.QuadPotential_Inv(inv), | ||
] | ||
if quadpotential.chol_available: | ||
pot = quadpotential.QuadPotential_Sparse(scipy.sparse.csc_matrix(cov)) | ||
pots.append(pot) | ||
for pot in pots: | ||
cov_ = np.cov(np.array([pot.random() for _ in range(10000)]).T) | ||
assert np.allclose(cov_, inv, atol=0.1) |