Datasets
Models
Diff of
/tests/edgepy/test_mglm.py
[000000]
..
[ea0fd6]
Switch to side-by-side view
--- a +++ b/tests/edgepy/test_mglm.py @@ -0,0 +1,290 @@ +import unittest + +import numpy as np + +from inmoose.edgepy import designAsFactor, mglmLevenberg, mglmOneGroup, mglmOneWay +from inmoose.utils import Factor, rnbinom + + +class test_mglm(unittest.TestCase): + def setUp(self): + y = np.array(rnbinom(80, size=5, mu=20, seed=42)).reshape((20, 4)) + self.y = np.vstack(([0, 0, 0, 0], [0, 0, 2, 2], y)) + self.group = np.array([1, 1, 2, 2]) + self.dispersion = 0.05 + self.lib_size = self.y.sum(axis=0) + + def test_mglmOneGroup(self): + j = self.group == 1 + y1 = mglmOneGroup( + self.y[:, j], dispersion=self.dispersion, offset=np.log(self.lib_size[j]) + ) + ref1 = np.array( + [ + -np.inf, + -np.inf, + -2.986383, + -3.485387, + -2.896575, + -3.676740, + -3.215694, + -3.056997, + -3.406943, + -3.423806, + -2.740691, + -3.226567, + -2.741902, + -2.648455, + -3.053310, + -2.971478, + -3.166231, + -2.754931, + -2.635907, + -2.494544, + -3.060501, + -3.192916, + ] + ) + self.assertTrue(np.allclose(y1, ref1, atol=1e-6, rtol=0)) + + j = self.group == 2 + y2 = mglmOneGroup( + self.y[:, j], dispersion=self.dispersion, offset=np.log(self.lib_size[j]) + ) + ref2 = np.array( + [ + -np.inf, + -5.420406, + -3.273714, + -3.506488, + -3.452560, + -2.762558, + -2.748622, + -2.780652, + -3.252326, + -2.851513, + -3.218678, + -3.277784, + -3.211742, + -3.407176, + -2.912961, + -2.475169, + -3.472935, + -3.057142, + -2.737961, + -2.636946, + -2.794320, + -3.094864, + ] + ) + self.assertTrue(np.allclose(y2, ref2, atol=1e-6, rtol=0)) + + # TODO also test with varying tolerance + def test_mglmLevenberg(self): + design = np.ones((self.y.shape[1], 1)) + coef_ref = np.array( + [ + np.nan, + 0, + 2.9704144655697013455, + 2.6026896854443828389, + 2.9704144655697013455, + 2.9831534913471302595, + 3.1463051320333654814, + 3.1986731175506815106, + 2.7725887222397807008, + 3.0081547935525478898, + 3.135494215929149675, + 2.8478121434773688847, + 3.135494215929149675, + 3.135494215929149675, + 3.1135153092103742267, + 3.4094961844768505443, + 2.788092908775746448, + 3.2188758248682005636, + 3.4094961844768505443, + 3.533686564708234279, + 3.1780538303479461959, + 2.9575110607337933288, + ] + ).reshape(self.y.shape[0], 1) + + fit_ref = np.array( + [ + [0.00, 0.00, 0.00, 0.00], + [1.00, 1.00, 1.00, 1.00], + [19.50, 19.50, 19.50, 19.50], + [13.50, 13.50, 13.50, 13.50], + [19.50, 19.50, 19.50, 19.50], + [19.75, 19.75, 19.75, 19.75], + [23.25, 23.25, 23.25, 23.25], + [24.50, 24.50, 24.50, 24.50], + [16.00, 16.00, 16.00, 16.00], + [20.25, 20.25, 20.25, 20.25], + [23.00, 23.00, 23.00, 23.00], + [17.25, 17.25, 17.25, 17.25], + [23.00, 23.00, 23.00, 23.00], + [23.00, 23.00, 23.00, 23.00], + [22.50, 22.50, 22.50, 22.50], + [30.25, 30.25, 30.25, 30.25], + [16.25, 16.25, 16.25, 16.25], + [25.00, 25.00, 25.00, 25.00], + [30.25, 30.25, 30.25, 30.25], + [34.25, 34.25, 34.25, 34.25], + [24.00, 24.00, 24.00, 24.00], + [19.25, 19.25, 19.25, 19.25], + ] + ) + + dev_ref = np.array( + [ + 0.0000000, + 5.5451767, + 6.4342458, + 3.6714251, + 40.1310924, + 19.0526072, + 7.5713327, + 22.3840805, + 0.8552751, + 21.3975973, + 9.7128441, + 2.8266627, + 21.4133481, + 11.9659419, + 3.8835243, + 8.5881018, + 1.6366558, + 49.5556034, + 17.5782768, + 4.3513363, + 4.2882649, + 1.0473488, + ] + ) + (coef, fit, dev, it, fail) = mglmLevenberg(self.y, design) + self.assertTrue(np.allclose(coef, coef_ref, atol=1e-15, rtol=0, equal_nan=True)) + self.assertTrue(np.allclose(fit, fit_ref, atol=1e-6, rtol=0)) + self.assertTrue(np.allclose(dev, dev_ref, atol=1e-6, rtol=0)) + self.assertEqual(it[0], 0) + self.assertEqual(np.unique(it[1:]), 1) + self.assertFalse(np.any(fail)) + + (coef, fit, dev, it, fail) = mglmLevenberg( + self.y, design, coef_start=np.ones((self.y.shape[0], design.shape[1])) + ) + self.assertTrue(np.allclose(coef, coef_ref, atol=1e-6, rtol=0, equal_nan=True)) + + with self.assertRaisesRegex(ValueError, expected_regex="no data"): + mglmLevenberg(np.ones(shape=(0, 1)), None) + with self.assertRaisesRegex( + ValueError, expected_regex="invalid start_method foo" + ): + mglmLevenberg(self.y, design, start_method="foo") + + +class test_mglmOneWay(unittest.TestCase): + def setUp(self): + y = np.array(rnbinom(80, size=5, mu=20, seed=42)).reshape((20, 4)) + self.y = np.vstack(([0, 0, 0, 0], [0, 0, 2, 2], y)) + self.group = np.array([1, 1, 2, 2]) + self.dispersion = 0.05 + self.lib_size = self.y.sum(axis=0) + + def test_designAsFactor(self): + f1 = designAsFactor(np.full((3, 4), [1, 2, 3, 4])) + f2 = designAsFactor(np.full((3, 4), [1, 2, 3, 4]).T) + self.assertTrue(np.array_equal(f1.__array__(), Factor([1, 1, 1]).__array__())) + self.assertTrue( + np.array_equal(f2.__array__(), Factor([1, 2, 3, 4]).__array__()) + ) + + def test_mglmOneWay(self): + coef_ref = np.array( + [ + -1.000000e08, + 0.000000e00, + 2.970414e00, + 2.602690e00, + 2.970414e00, + 2.983153e00, + 3.146305e00, + 3.198673e00, + 2.772589e00, + 3.008155e00, + 3.135494e00, + 2.847812e00, + 3.135494e00, + 3.135494e00, + 3.113515e00, + 3.409496e00, + 2.788093e00, + 3.218876e00, + 3.409496e00, + 3.533687e00, + 3.178054e00, + 2.957511e00, + ] + ).reshape(self.y.shape[0], 1) + (coef, fit) = mglmOneWay(self.y) + self.assertTrue(np.allclose(coef, coef_ref, atol=1e-6, rtol=0)) + + design = np.array([[1, 0], [1, 0], [0, 1], [0, 1]]) + coef_ref = np.array( + [ + [-1.000000e08, -1.000000e08], + [-1.000000e08, 6.931472e-01], + [3.091042e00, 2.833213e00], + [2.602690e00, 2.602690e00], + [3.198673e00, 2.674149e00], + [2.397895e00, 3.349904e00], + [2.862201e00, 3.367296e00], + [3.044522e00, 3.332205e00], + [2.674149e00, 2.862201e00], + [2.674149e00, 3.258097e00], + [3.332205e00, 2.890372e00], + [2.862201e00, 2.833213e00], + [3.332205e00, 2.890372e00], + [3.433987e00, 2.708050e00], + [3.020425e00, 3.198673e00], + [3.113515e00, 3.637586e00], + [2.917771e00, 2.639057e00], + [3.349904e00, 3.068053e00], + [3.433987e00, 3.384390e00], + [3.583519e00, 3.481240e00], + [3.020425e00, 3.314186e00], + [2.890372e00, 3.020425e00], + ] + ) + (coef, fit) = mglmOneWay(self.y, design=design) + self.assertTrue(np.allclose(coef, coef_ref, atol=1e-6, rtol=0)) + + design = np.array([[2, 0], [2, 0], [0, 2], [0, 2]]) + coef_ref = np.array( + [ + [-5.000000e07, -5.000000e07], + [-5.000000e07, 3.465736e-01], + [1.545521e00, 1.416607e00], + [1.301345e00, 1.301345e00], + [1.599337e00, 1.337074e00], + [1.198948e00, 1.674952e00], + [1.431100e00, 1.683648e00], + [1.522261e00, 1.666102e00], + [1.337074e00, 1.431100e00], + [1.337074e00, 1.629048e00], + [1.666102e00, 1.445186e00], + [1.431100e00, 1.416607e00], + [1.666102e00, 1.445186e00], + [1.716994e00, 1.354025e00], + [1.510212e00, 1.599337e00], + [1.556758e00, 1.818793e00], + [1.458885e00, 1.319529e00], + [1.674952e00, 1.534026e00], + [1.716994e00, 1.692195e00], + [1.791759e00, 1.740620e00], + [1.510212e00, 1.657093e00], + [1.445186e00, 1.510212e00], + ] + ) + (coef, fit) = mglmOneWay(self.y, design=design) + self.assertTrue(np.allclose(coef, coef_ref, atol=1e-6, rtol=0))