def test_cartesian():
# Check if cartesian product delivers the right results
axes = (np.array([1, 2, 3]), np.array([4, 5]), np.array([6, 7]))
true_out = np.array([[1, 4, 6],
[1, 4, 7],
[1, 5, 6],
[1, 5, 7],
[2, 4, 6],
[2, 4, 7],
[2, 5, 6],
[2, 5, 7],
[3, 4, 6],
[3, 4, 7],
[3, 5, 6],
[3, 5, 7]])
out = cartesian(axes)
assert_array_equal(true_out, out)
# check single axis
x = np.arange(3)
assert_array_equal(x[:, np.newaxis], cartesian((x,)))
python类cartesian()的实例源码
def test_cgauss_likelihood():
mu = np.array([0], dtype='float')
sigma = np.array([2], dtype='float')
x = np.linspace(-1, 2, 2)
lapse = np.array([0], dtype='float')
parameters = cartesian((mu, sigma, lapse, x))
proportionMethod = PsiMarginal.pf(parameters, psyfun='cGauss')
samples = np.random.normal(mu, sigma, (200000, 1))
proportionSamples = np.empty([2, ])
proportionSamples[0] = np.mean(samples <= x[0]) # cdf is p(X<=x), compute this through sampling to check likelihood
proportionSamples[1] = np.mean(samples <= x[1])
np.testing.assert_almost_equal(proportionSamples, proportionMethod, decimal=2) == 1
def compute_overlap(mat1, mat2):
s1 = mat1.shape[0]
s2 = mat2.shape[0]
area1 = (mat1[:, 2] - mat1[:, 0]) * (mat1[:, 3] - mat1[:, 1])
if mat2.shape[1] == 5:
area2 = mat2[:, 4]
else:
area2 = (mat2[:, 2] - mat2[:, 0]) * (mat2[:, 3] - mat2[:, 1])
x1 = cartesian([mat1[:, 0], mat2[:, 0]])
x1 = np.amax(x1, axis=1)
x2 = cartesian([mat1[:, 2], mat2[:, 2]])
x2 = np.amin(x2, axis=1)
com_zero = np.zeros(x2.shape[0])
w = x2 - x1
w = w - 1
w = np.maximum(com_zero, w)
y1 = cartesian([mat1[:, 1], mat2[:, 1]])
y1 = np.amax(y1, axis=1)
y2 = cartesian([mat1[:, 3], mat2[:, 3]])
y2 = np.amin(y2, axis=1)
h = y2 - y1
h = h - 1
h = np.maximum(com_zero, h)
oo = w * h
aa = cartesian([area1[:], area2[:]])
aa = np.sum(aa, axis=1)
ooo = oo / (aa - oo)
overlap = np.transpose(ooo.reshape(s1, s2), (1, 0))
return overlap
def weights(dim, degree):
# 1D sigma-points (x) and weights (w)
x, w = hermegauss(degree)
# hermegauss() provides weights that cause posdef errors
w = factorial(degree) / (degree ** 2 * hermeval(x, [0] * (degree - 1) + [1]) ** 2)
return np.prod(cartesian([w] * dim), axis=1)
def unit_sigma_points(dim, degree):
# 1D sigma-points (x) and weights (w)
x, w = hermegauss(degree)
# nD sigma-points by cartesian product
return cartesian([x] * dim).T # column/sigma-point
