From e1bfc4b48235792ab1ca4e128bd0f6ad61d799fe Mon Sep 17 00:00:00 2001 From: Chris Tralie Date: Tue, 14 Jul 2020 21:27:37 -0400 Subject: [PATCH 1/3] Adding test from #14 --- cechmate/solver.py | 2 +- test/test_alpha.py | 38 ++++++++++++++++++++++++++++++++++++++ 2 files changed, 39 insertions(+), 1 deletion(-) create mode 100644 test/test_alpha.py diff --git a/cechmate/solver.py b/cechmate/solver.py index 7e9bf14..019590e 100644 --- a/cechmate/solver.py +++ b/cechmate/solver.py @@ -87,7 +87,7 @@ def _simplices_to_sparse_pivot_column(ordered_simplices, verbose=False): fidxs = np.array(list(fidxs)) fidxs = tuple(idxs[fidxs]) if not fidxs in idxs2order: - print( + raise Exception( "Error: Not a proper filtration: %s added before %s" % (idxs, fidxs) ) diff --git a/test/test_alpha.py b/test/test_alpha.py new file mode 100644 index 0000000..57f4329 --- /dev/null +++ b/test/test_alpha.py @@ -0,0 +1,38 @@ +import pytest + +import numpy as np +from cechmate import Alpha + +@pytest.fixture +def triangle(): + x = np.array([ + [0, 0.0], + [1, 1.0], + [0, 1.0], + ]) + + return x + +def test_triangle(triangle): + """ Expect 3 vertices, 3 edges, and a triangle + + """ + a = Alpha(2).build(triangle) + + assert len(a) == 7 + + vertices = [s for s in a if len(s[0]) == 1] + edges = [s for s in a if len(s[0]) == 2] + triangles = [s for s in a if len(s[0]) == 3] + + assert len(vertices) == 3 + assert len(edges) == 3 + assert len(triangles) == 1 + +def test_precision(): + X = np.array([4.148168134442850770e-16,2.579509799999999853e+00,3.403597400000000217e+00,1.846206783921459376e-15,1.148050980000000010e+01,1.237459740000000075e+01,6.751000000000002110e+00,1.148050980000000010e+01,3.403597400000001105e+00,6.751000000000000334e+00,2.579509799999999853e+00,1.237459740000000075e+01,2.447963127510063373e-15,1.522249019999999931e+01,3.403597400000001105e+00,1.016573157032889123e-15,6.321490199999999504e+00,1.237459740000000075e+01,6.751000000000001222e+00,6.321490199999999504e+00,3.403597400000001105e+00,6.751000000000002998e+00,1.522249019999999931e+01,1.237459740000000075e+01,2.447963127510063373e-15,1.522249019999999931e+01,1.453840260000000129e+01,1.016573157032889123e-15,6.321490199999999504e+00,5.567402600000001200e+00,6.751000000000001222e+00,6.321490199999999504e+00,1.453840260000000129e+01,6.751000000000002998e+00,1.522249019999999931e+01,5.567402600000002089e+00,4.148168134442850770e-16,2.579509799999999853e+00,1.453840260000000129e+01,1.846206783921459376e-15,1.148050980000000010e+01,5.567402600000001200e+00,6.751000000000002110e+00,1.148050980000000010e+01,1.453840260000000129e+01,6.751000000000000334e+00,2.579509799999999853e+00,5.567402600000001200e+00,1.806567600000000384e+00,1.313787599999999944e+00,4.837163199999999996e+00,1.806567600000001717e+00,1.021478759999999930e+01,1.380816320000000275e+01,8.557567600000002273e+00,1.021478759999999930e+01,4.837163200000000884e+00,8.557567600000000496e+00,1.313787599999999944e+00,1.380816320000000097e+01,1.169543239999999962e+01,1.313787599999999944e+00,4.837163200000000884e+00,1.169543240000000139e+01,1.021478759999999930e+01,1.380816320000000275e+01,4.944432400000001948e+00,1.021478759999999930e+01,4.837163200000000884e+00,4.944432400000000172e+00,1.313787599999999944e+00,1.380816320000000097e+01,1.806567600000002827e+00,1.648821239999999833e+01,4.837163200000000884e+00,1.806567600000001272e+00,7.587212400000000301e+00,1.380816320000000097e+01,8.557567600000002273e+00,7.587212400000000301e+00,4.837163200000000884e+00,8.557567600000002273e+00,1.648821239999999833e+01,1.380816320000000275e+01,1.169543240000000139e+01,1.648821239999999833e+01,4.837163200000001773e+00,1.169543240000000139e+01,7.587212400000000301e+00,1.380816320000000275e+01,4.944432400000001060e+00,7.587212400000000301e+00,4.837163200000000884e+00,4.944432400000002836e+00,1.648821239999999833e+01,1.380816320000000275e+01,1.169543240000000139e+01,1.648821239999999833e+01,1.310483680000000106e+01,1.169543240000000139e+01,7.587212400000000301e+00,4.133836800000000977e+00,4.944432400000001060e+00,7.587212400000000301e+00,1.310483679999999929e+01,4.944432400000002836e+00,1.648821239999999833e+01,4.133836800000000977e+00,1.806567600000002827e+00,1.648821239999999833e+01,1.310483679999999929e+01,1.806567600000001272e+00,7.587212400000000301e+00,4.133836800000000089e+00,8.557567600000002273e+00,7.587212400000000301e+00,1.310483679999999929e+01,8.557567600000002273e+00,1.648821239999999833e+01,4.133836800000001865e+00,1.169543239999999962e+01,1.313787599999999944e+00,1.310483679999999929e+01,1.169543240000000139e+01,1.021478759999999930e+01,4.133836800000000977e+00,4.944432400000001948e+00,1.021478759999999930e+01,1.310483679999999929e+01,4.944432400000000172e+00,1.313787599999999944e+00,4.133836800000000089e+00,1.806567600000000384e+00,1.313787599999999944e+00,1.310483679999999929e+01,1.806567600000001717e+00,1.021478759999999930e+01,4.133836800000000977e+00,8.557567600000002273e+00,1.021478759999999930e+01,1.310483679999999929e+01,8.557567600000000496e+00,1.313787599999999944e+00,4.133836800000000089e+00,2.446562400000000359e+00,3.521235599999999799e+00,3.556104400000000165e+00,2.446562400000001691e+00,1.242223559999999871e+01,1.252710439999999892e+01,9.197562400000002469e+00,1.242223559999999871e+01,3.556104400000001053e+00,9.197562400000000693e+00,3.521235599999999799e+00,1.252710439999999892e+01,1.105543759999999942e+01,3.521235599999999799e+00,3.556104400000000609e+00,1.105543760000000120e+01,1.242223559999999871e+01,1.252710440000000069e+01,4.304437600000001751e+00,1.242223559999999871e+01,3.556104400000000609e+00,4.304437600000000863e+00,3.521235599999999799e+00,1.252710439999999892e+01,2.446562400000002135e+00,1.428076440000000069e+01,3.556104400000000609e+00,2.446562400000000803e+00,5.379764400000000002e+00,1.252710439999999892e+01,9.197562400000000693e+00,5.379764400000000002e+00,3.556104400000000609e+00,9.197562400000002469e+00,1.428076440000000069e+01,1.252710440000000069e+01,1.105543760000000120e+01,1.428076440000000069e+01,3.556104400000001498e+00,1.105543759999999942e+01,5.379764400000000002e+00,1.252710440000000069e+01,4.304437600000000863e+00,5.379764400000000002e+00,3.556104400000000609e+00,4.304437600000002639e+00,1.428076440000000069e+01,1.252710440000000069e+01,1.105543760000000120e+01,1.428076440000000069e+01,1.438589560000000311e+01,1.105543759999999942e+01,5.379764400000000002e+00,5.414895600000001252e+00,4.304437600000000863e+00,5.379764400000000002e+00,1.438589560000000134e+01,4.304437600000002639e+00,1.428076440000000069e+01,5.414895600000001252e+00,2.446562400000002135e+00,1.428076440000000069e+01,1.438589560000000311e+01,2.446562400000000803e+00,5.379764400000000002e+00,5.414895600000000364e+00,9.197562400000000693e+00,5.379764400000000002e+00,1.438589560000000311e+01,9.197562400000002469e+00,1.428076440000000069e+01,5.414895600000002140e+00,1.105543759999999942e+01,3.521235599999999799e+00,1.438589560000000311e+01,1.105543760000000120e+01,1.242223559999999871e+01,5.414895600000002140e+00,4.304437600000001751e+00,1.242223559999999871e+01,1.438589560000000311e+01,4.304437600000000863e+00,3.521235599999999799e+00,5.414895600000000364e+00,2.446562400000000359e+00,3.521235599999999799e+00,1.438589560000000134e+01,2.446562400000001691e+00,1.242223559999999871e+01,5.414895600000001252e+00,9.197562400000002469e+00,1.242223559999999871e+01,1.438589560000000311e+01,9.197562400000000693e+00,3.521235599999999799e+00,5.414895600000001252e+00,1.972642200000000345e+00,1.313787599999999944e+00,2.212248600000000565e+00,1.972642200000001678e+00,1.021478759999999930e+01,1.118324860000000065e+01,8.723642200000002234e+00,1.021478759999999930e+01,2.212248600000001453e+00,8.723642200000000457e+00,1.313787599999999944e+00,1.118324860000000065e+01,1.152935779999999966e+01,1.313787599999999944e+00,2.212248600000001009e+00,1.152935780000000143e+01,1.021478759999999930e+01,1.118324860000000065e+01,4.778357800000001987e+00,1.021478759999999930e+01,2.212248600000001009e+00,4.778357800000000211e+00,1.313787599999999944e+00,1.118324860000000065e+01,1.972642200000002788e+00,1.648821239999999833e+01,2.212248600000001453e+00,1.972642200000001234e+00,7.587212400000000301e+00,1.118324860000000065e+01,8.723642200000002234e+00,7.587212400000000301e+00,2.212248600000001009e+00,8.723642200000002234e+00,1.648821239999999833e+01,1.118324860000000065e+01,1.152935780000000143e+01,1.648821239999999833e+01,2.212248600000001897e+00,1.152935780000000143e+01,7.587212400000000301e+00,1.118324860000000065e+01,4.778357800000001099e+00,7.587212400000000301e+00,2.212248600000001009e+00,4.778357800000002875e+00,1.648821239999999833e+01,1.118324860000000065e+01,1.152935780000000143e+01,1.648821239999999833e+01,1.572975140000000316e+01,1.152935780000000143e+01,7.587212400000000301e+00,6.758751400000001297e+00,4.778357800000001099e+00,7.587212400000000301e+00,1.572975140000000138e+01,4.778357800000002875e+00,1.648821239999999833e+01,6.758751400000001297e+00,1.972642200000002788e+00,1.648821239999999833e+01,1.572975140000000138e+01,1.972642200000001234e+00,7.587212400000000301e+00,6.758751400000000409e+00,8.723642200000002234e+00,7.587212400000000301e+00,1.572975140000000138e+01,8.723642200000002234e+00,1.648821239999999833e+01,6.758751400000001297e+00,1.152935779999999966e+01,1.313787599999999944e+00,1.572975140000000138e+01,1.152935780000000143e+01,1.021478759999999930e+01,6.758751400000001297e+00,4.778357800000001987e+00,1.021478759999999930e+01,1.572975140000000138e+01,4.778357800000000211e+00,1.313787599999999944e+00,6.758751400000000409e+00,1.972642200000000345e+00,1.313787599999999944e+00,1.572975140000000138e+01,1.972642200000001678e+00,1.021478759999999930e+01,6.758751400000000409e+00,8.723642200000002234e+00,1.021478759999999930e+01,1.572975140000000138e+01,8.723642200000000457e+00,1.313787599999999944e+00,6.758751400000000409e+00,4.236927600000000460e+00,0.000000000000000000e+00,1.971825800000000184e+00,4.236927600000002236e+00,8.900999999999999801e+00,1.094282580000000138e+01,1.098792760000000257e+01,8.900999999999999801e+00,1.971825800000001294e+00,1.098792760000000079e+01,0.000000000000000000e+00,1.094282580000000138e+01,9.265072399999999320e+00,0.000000000000000000e+00,1.971825800000000628e+00,9.265072400000001096e+00,8.900999999999999801e+00,1.094282580000000138e+01,2.514072400000001206e+00,8.900999999999999801e+00,1.971825800000000628e+00,2.514072399999999874e+00,0.000000000000000000e+00,1.094282579999999960e+01,9.265072399999999320e+00,0.000000000000000000e+00,1.597017420000000065e+01,9.265072400000001096e+00,8.900999999999999801e+00,6.999174200000001456e+00,2.514072400000001206e+00,8.900999999999999801e+00,1.597017420000000065e+01,2.514072399999999874e+00,0.000000000000000000e+00,6.999174200000000567e+00,4.236927600000000460e+00,0.000000000000000000e+00,1.597017420000000065e+01,4.236927600000002236e+00,8.900999999999999801e+00,6.999174200000000567e+00,1.098792760000000257e+01,8.900999999999999801e+00,1.597017420000000065e+01,1.098792760000000079e+01,0.000000000000000000e+00,6.999174200000000567e+00,2.492469199999999940e+00,0.000000000000000000e+00,1.526197213876682126e-16,2.492469200000001273e+00,8.900999999999999801e+00,8.971000000000000085e+00,9.243469200000001607e+00,8.900999999999999801e+00,1.111028306400386806e-15,9.243469199999999830e+00,0.000000000000000000e+00,8.971000000000000085e+00,1.100953080000000028e+01,0.000000000000000000e+00,6.741393327167100174e-16,1.100953080000000206e+01,8.900999999999999801e+00,8.971000000000001862e+00,4.258530800000002614e+00,8.900999999999999801e+00,8.057888636250504791e-16,4.258530800000000838e+00,0.000000000000000000e+00,8.971000000000000085e+00,0.000000000000000000e+00,0.000000000000000000e+00,6.224079800000000162e+00,1.431389970477174250e-15,8.900999999999999801e+00,1.519507980000000025e+01,6.751000000000002110e+00,8.900999999999999801e+00,6.224079800000001050e+00,6.751000000000000334e+00,0.000000000000000000e+00,1.519507980000000025e+01,0.000000000000000000e+00,0.000000000000000000e+00,1.171792020000000001e+01,1.431389970477174250e-15,8.900999999999999801e+00,2.746920200000000811e+00,6.751000000000002110e+00,8.900999999999999801e+00,1.171792020000000178e+01,6.751000000000000334e+00,0.000000000000000000e+00,2.746920200000000811e+00,1.555430400000000546e+00,2.180744999999999933e+00,3.504072600000000204e+00,1.555430400000001878e+00,1.108174500000000151e+01,1.247507260000000073e+01,8.306430400000001768e+00,1.108174500000000151e+01,3.504072600000001092e+00,8.306430399999999992e+00,2.180744999999999933e+00,1.247507260000000073e+01,1.194656960000000012e+01,2.180744999999999933e+00,3.504072600000001092e+00,1.194656960000000190e+01,1.108174500000000151e+01,1.247507260000000251e+01,5.195569600000002453e+00,1.108174500000000151e+01,3.504072600000001092e+00,5.195569600000000676e+00,2.180744999999999933e+00,1.247507260000000073e+01,1.555430400000002544e+00,1.562125499999999789e+01,3.504072600000001092e+00,1.555430400000001212e+00,6.720254999999999868e+00,1.247507260000000073e+01,8.306430400000001768e+00,6.720254999999999868e+00,3.504072600000001092e+00,8.306430400000001768e+00,1.562125499999999789e+01,1.247507260000000251e+01,1.194656960000000190e+01,1.562125499999999789e+01,3.504072600000001536e+00,1.194656960000000190e+01,6.720254999999999868e+00,1.247507260000000251e+01,5.195569600000001564e+00,6.720254999999999868e+00,3.504072600000000648e+00,5.195569600000003341e+00,1.562125499999999789e+01,1.247507260000000251e+01,1.194656960000000190e+01,1.562125499999999789e+01,1.443792740000000130e+01,1.194656960000000190e+01,6.720254999999999868e+00,5.466927400000000326e+00,5.195569600000001564e+00,6.720254999999999868e+00,1.443792739999999952e+01,5.195569600000003341e+00,1.562125499999999789e+01,5.466927400000001214e+00,1.555430400000002544e+00,1.562125499999999789e+01,1.443792740000000130e+01,1.555430400000001212e+00,6.720254999999999868e+00,5.466927400000000326e+00,8.306430400000001768e+00,6.720254999999999868e+00,1.443792740000000130e+01,8.306430400000001768e+00,1.562125499999999789e+01,5.466927400000001214e+00,1.194656960000000012e+01,2.180744999999999933e+00,1.443792740000000130e+01,1.194656960000000190e+01,1.108174500000000151e+01,5.466927400000001214e+00,5.195569600000002453e+00,1.108174500000000151e+01,1.443792740000000130e+01,5.195569600000000676e+00,2.180744999999999933e+00,5.466927400000000326e+00,1.555430400000000546e+00,2.180744999999999933e+00,1.443792739999999952e+01,1.555430400000001878e+00,1.108174500000000151e+01,5.466927400000000326e+00,8.306430400000001768e+00,1.108174500000000151e+01,1.443792740000000130e+01,8.306430399999999992e+00,2.180744999999999933e+00,5.466927400000000326e+00,2.670695600000000169e+00,0.000000000000000000e+00,1.596838000000000202e+00,2.670695600000001502e+00,8.900999999999999801e+00,1.056783800000000006e+01,9.421695600000001392e+00,8.900999999999999801e+00,1.596838000000001090e+00,9.421695599999999615e+00,0.000000000000000000e+00,1.056783800000000006e+01,1.083130440000000050e+01,0.000000000000000000e+00,1.596838000000000646e+00,1.083130440000000227e+01,8.900999999999999801e+00,1.056783800000000006e+01,4.080304400000001941e+00,8.900999999999999801e+00,1.596838000000000646e+00,4.080304400000000165e+00,0.000000000000000000e+00,1.056783800000000006e+01,1.083130440000000050e+01,0.000000000000000000e+00,1.634516200000000197e+01,1.083130440000000227e+01,8.900999999999999801e+00,7.374162000000001882e+00,4.080304400000001941e+00,8.900999999999999801e+00,1.634516200000000197e+01,4.080304400000000165e+00,0.000000000000000000e+00,7.374162000000000994e+00,2.670695600000000169e+00,0.000000000000000000e+00,1.634516200000000197e+01,2.670695600000001502e+00,8.900999999999999801e+00,7.374162000000000994e+00,9.421695600000001392e+00,8.900999999999999801e+00,1.634516200000000197e+01,9.421695599999999615e+00,0.000000000000000000e+00,7.374162000000000994e+00,1.531126800000000010e+00,0.000000000000000000e+00,5.723498000000000197e+00,1.531126800000001342e+00,8.900999999999999801e+00,1.469449799999999939e+01,8.282126800000002120e+00,8.900999999999999801e+00,5.723498000000001085e+00,8.282126800000000344e+00,0.000000000000000000e+00,1.469449799999999939e+01,1.197087320000000155e+01,0.000000000000000000e+00,5.723498000000001085e+00,1.197087320000000332e+01,8.900999999999999801e+00,1.469449800000000117e+01,5.219873200000002100e+00,8.900999999999999801e+00,5.723498000000001085e+00,5.219873200000000324e+00,0.000000000000000000e+00,1.469449799999999939e+01,1.197087320000000155e+01,0.000000000000000000e+00,1.221850200000000086e+01,1.197087320000000332e+01,8.900999999999999801e+00,3.247502000000001221e+00,5.219873200000002100e+00,8.900999999999999801e+00,1.221850200000000264e+01,5.219873200000000324e+00,0.000000000000000000e+00,3.247502000000000333e+00,1.531126800000000010e+00,0.000000000000000000e+00,1.221850200000000086e+01,1.531126800000001342e+00,8.900999999999999801e+00,3.247502000000000777e+00,8.282126800000002120e+00,8.900999999999999801e+00,1.221850200000000264e+01,8.282126800000000344e+00,0.000000000000000000e+00,3.247502000000000333e+00,3.375500000000001055e+00,4.450499999999999901e+00,4.485500000000000931e+00,3.375500000000002387e+00,1.335149999999999970e+01,1.345650000000000190e+01,1.012650000000000183e+01,1.335149999999999970e+01,4.485500000000001819e+00,1.012650000000000006e+01,4.450499999999999901e+00,1.345650000000000190e+01,1.012650000000000006e+01,4.450499999999999901e+00,4.485500000000000931e+00,1.012650000000000183e+01,1.335149999999999970e+01,1.345650000000000190e+01,3.375500000000002387e+00,1.335149999999999970e+01,4.485500000000000931e+00,3.375500000000001055e+00,4.450499999999999901e+00,1.345650000000000013e+01]) + X = np.reshape(X, (int(X.size/3), 3)) + alpha = Alpha() + alpha_filtration = alpha.build(X) + dgms = alpha.diagrams(alpha_filtration) + assert(len(dgms) == 4) \ No newline at end of file From db2c3e4089aa768d0e01f39911621a086780ada5 Mon Sep 17 00:00:00 2001 From: Chris Tralie Date: Tue, 14 Jul 2020 22:05:26 -0400 Subject: [PATCH 2/3] Improving numerical precision and fixing simplex condition bug --- cechmate/filtrations/alpha.py | 66 +++++++++++++++++++++-------------- 1 file changed, 40 insertions(+), 26 deletions(-) diff --git a/cechmate/filtrations/alpha.py b/cechmate/filtrations/alpha.py index 8b33fc0..d878ab8 100644 --- a/cechmate/filtrations/alpha.py +++ b/cechmate/filtrations/alpha.py @@ -2,7 +2,7 @@ import time import numpy as np -import numpy.linalg as linalg +from numpy import linalg from scipy import spatial from .base import BaseFiltration @@ -25,6 +25,7 @@ class Alpha(BaseFiltration): >>> diagrams = r.diagrams(simplices) """ + MIN_DET = 1e-10 def build(self, X): """ @@ -61,35 +62,40 @@ def build(self, X): tic = time.time() filtration = {} - simplices_bydim = {} for dim in range(maxdim + 2, 1, -1): - simplices_bydim[dim] = [] for s in range(delaunay_faces.shape[0]): simplex = delaunay_faces[s, :] for sigma in itertools.combinations(simplex, dim): sigma = tuple(sorted(sigma)) - simplices_bydim[dim].append(sigma) if not sigma in filtration: - filtration[sigma] = self._get_circumcenter(X[sigma, :])[1] - for i in range(dim): - # Propagate alpha filtration value - tau = sigma[0:i] + sigma[i + 1 : :] - if tau in filtration: - filtration[tau] = min(filtration[tau], filtration[sigma]) - elif len(tau) > 1: - # If Tau is not empty - xtau, rtauSqr = self._get_circumcenter(X[tau, :]) - if np.sum((X[sigma[i], :] - xtau) ** 2) < rtauSqr: - filtration[tau] = filtration[sigma] - for f in filtration: - filtration[f] = np.sqrt(filtration[f]) + rSqr = self._get_circumcenter(X[sigma, :])[1] + if np.isfinite(rSqr): + filtration[sigma] = rSqr + if sigma in filtration: + for i in range(dim): # Propagate alpha filtration value + tau = sigma[0:i] + sigma[i+1::] + if tau in filtration: + filtration[tau] = min(filtration[tau], filtration[sigma]) + elif len(tau) > 1 and sigma in filtration: + # If Tau is not empty + xtau, rtauSqr = self._get_circumcenter(X[tau, :]) + if np.sum((X[sigma[i], :] - xtau) ** 2) < rtauSqr: + filtration[tau] = filtration[sigma] + # Convert from squared radii to radii + for sigma in filtration: + filtration[sigma] = np.sqrt(filtration[sigma]) ## Step 2: Take care of numerical artifacts that may result ## in simplices with greater filtration values than their co-faces - for dim in range(maxdim + 2, 2, -1): - for sigma in simplices_bydim[dim]: - for i in range(dim): - tau = sigma[0:i] + sigma[i + 1 : :] + simplices_bydim = [set([]) for i in range(maxdim+2)] + for simplex in filtration.keys(): + simplices_bydim[len(simplex)-1].add(simplex) + simplices_bydim = simplices_bydim[2::] + simplices_bydim.reverse() + for simplices_dim in simplices_bydim: + for sigma in simplices_dim: + for i in range(len(sigma)): + tau = sigma[0:i] + sigma[i+1::] if filtration[tau] > filtration[sigma]: filtration[tau] = filtration[sigma] @@ -127,6 +133,7 @@ def _get_circumcenter(self, X): (SC3) If there are more points than the ambient dimension plus one it returns (np.nan, np.nan) """ + X0 = np.array(X) if X.shape[0] == 2: # Special case of an edge, which is very simple dX = X[1, :] - X[0, :] @@ -142,8 +149,7 @@ def _get_circumcenter(self, X): # Transform arrays for PCA for SC1 (points in higher ambient dimension) muV = np.array([]) V = np.array([]) - if X.shape[0] < X.shape[1] + 1: - # SC1: Do PCA down to NPoints-1 + if X.shape[0] < X.shape[1] + 1: # SC1: Do PCA down to NPoints-1 muV = np.mean(X, 0) XCenter = X - muV _, V = linalg.eigh((XCenter.T).dot(XCenter)) @@ -151,19 +157,27 @@ def _get_circumcenter(self, X): X = XCenter.dot(V) muX = np.mean(X, 0) D = np.ones((X.shape[0], X.shape[0] + 1)) - # Subtract off centroid for numerical stability - D[:, 1:-1] = X - muX + # Subtract off centroid and scale down for numerical stability + Y = X - muX + scaleSqr = np.max(np.sum(Y**2, 1)) + scaleSqr = 1 + scale = np.sqrt(scaleSqr) + Y /=scale + + D[:, 1:-1] = Y D[:, 0] = np.sum(D[:, 1:-1] ** 2, 1) minor = lambda A, j: A[ :, np.concatenate((np.arange(j), np.arange(j + 1, A.shape[1]))) ] dxs = np.array([linalg.det(minor(D, i)) for i in range(1, D.shape[1] - 1)]) alpha = linalg.det(minor(D, 0)) - if np.abs(alpha) > 0: + if np.abs(alpha) > Alpha.MIN_DET: signs = (-1) ** np.arange(len(dxs)) x = dxs * signs / (2 * alpha) + muX # Add back centroid gamma = ((-1) ** len(dxs)) * linalg.det(minor(D, D.shape[1] - 1)) rSqr = (np.sum(dxs ** 2) + 4 * alpha * gamma) / (4 * alpha * alpha) + x *= scale + rSqr *= scaleSqr if V.size > 0: # Transform back to ambient if SC1 x = x.dot(V.T) + muV From 637dd720ede3f2b8dc030dc21efcc0beca8f4473 Mon Sep 17 00:00:00 2001 From: Chris Tralie Date: Tue, 14 Jul 2020 22:09:47 -0400 Subject: [PATCH 3/3] Mistake in test, updating version, fixing readme link --- README.md | 2 +- cechmate/_version.py | 2 +- test/test_alpha.py | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/README.md b/README.md index 8ae5470..ecea633 100644 --- a/README.md +++ b/README.md @@ -37,4 +37,4 @@ To contribute please fork the project make your changes and submit a pull reques ## Documentation -Check out complete documentation at [cechmate.scikit-tda.org](http://www.cechmate.scikit-tda.org/) \ No newline at end of file +Check out complete documentation at [cechmate.scikit-tda.org](https://cechmate.scikit-tda.org/) \ No newline at end of file diff --git a/cechmate/_version.py b/cechmate/_version.py index a73339b..00ec2dc 100644 --- a/cechmate/_version.py +++ b/cechmate/_version.py @@ -1 +1 @@ -__version__ = "0.0.8" +__version__ = "0.0.9" diff --git a/test/test_alpha.py b/test/test_alpha.py index 57f4329..8485ac1 100644 --- a/test/test_alpha.py +++ b/test/test_alpha.py @@ -35,4 +35,4 @@ def test_precision(): alpha = Alpha() alpha_filtration = alpha.build(X) dgms = alpha.diagrams(alpha_filtration) - assert(len(dgms) == 4) \ No newline at end of file + assert(len(dgms) == 3)